We have succeeded in development of a simulation specialized for GI POF. In this study, we investigated the propagation characteristics of GI POF by use of this simulation. Propagation properties of multi-mode optical fibers can be calculated by the scalar-wave equation derived from Maxwell's equations. However, calculated impulse response disagrees with measured results. The factors of this disagreement have been generally explained as mode coupling and differential mode attenuation. These effects can be calculated by the power flow equation, as it has been applied for analysis of glass optical fibers and step-index polymer optical fibers. In this study, we applied the power flow equation to the graded-index polymer optical fiber (GI POF). The equation contains several parameters: propagation constants, coupling coefficients, and attenuation coefficient. In order to define these parameters, we fabricated poly methyl methacrylate (PMMA) based GI POF. Propagation constants of the GI POF were calculated by use of the finite-element method. Coupling and attenuation coefficients were estimated based on comparisons of measurements with simulation of differential mode attenuation and differential mode delay. We assigned these values to the power flow equation and solved it by use of the finite difference method. As a result, bandwidth characteristics calculated by this simulation well agreed with measurements. Moreover, it was found that the effect of mode coupling on impulse response of GI POF was more influential than that of differential mode attenuation and that higher modes were subject to mode coupling than lower modes and they were coupled into lower mode.