The solution of initial value problems provide the state history for a given dynamic system, for prescribed initial conditions. Existing methods for solving these problems have not been very successful in exploiting parallel computation architectures, mainly because most of the integration methods implemented on parallel machines are only modified versions of forward integration approaches, which are typically poorly suited for parallel computation. This paper propose parallel-structured modified Chebyshev-Picard iteration (MCPI) methods, which iteratively refine estimates of the solutions until the iteration converges. Using Chebyshev polynomials as the orthogonal approximation basis, it is straightforward to distribute the computation of force functions and polynomial coefficients to different processors. A vector-matrix form is introduced that makes the methods computationally efficient. The power of the methods is illustrated through satellite motion propagation problems. Compared with a Runge-Kutta 4-5 forward integration method implemented in MATLAB, the proposed methods generate solutions with improved accuracy as well as several orders of magnitude speedup, even prior to parallel implementation. Allowing to only integrate position states or perturbation motion achieve further speedup. Parallel realization of the methods is implemented using a graphics processing unit to provide an inexpensive parallel computation architecture. Significant further speedup is achieved from the parallel implementation.