Abstract This paper extends previous work on parallel-structured Modified Chebyshev Picard Iteration (MCPI) Methods. The MCPI approach iteratively refines path approximation of the state trajectory for smooth nonlinear dynamical systems and this paper shows that the approach is especially suitable for initial value problems of astrodynamics. Using Chebyshev polynomials, as the orthogonal approximation basis, it is straightforward to distribute the computation of force functions needed in MCPI to generate the polynomial coefficients (approximating the path iterations) to different processors. Combining Chebyshev polynomials with Picard iteration, MCPI methods iteratively refines path estimates over large time intervals chosen to be within the domain of convergence of Picard iteration. The developed vector-matrix form makes MCPI methods computationally efficient and a more systematic approach is given, leading to a modest correction to results in the published dissertation by Bai. The power of MCPI methods for solving IVPs is clearly illustrated using a simple nonlinear differential equation with a known analytical solution. Compared with the most common integration scheme, the standard Runge-Kutta 4-5 method as implemented in MATLAB, MCPI methods generate solutions with better accuracy as well as orders of magnitude speedups, on a serial machine. MCPI performance is also compared to state of the art integrators such as the Runge-Kutta Nystrom 12(10) methods applied to the relevant orbit mechanics problems. The MCPI method is shown to be well-suited to solving these problems in serial processors with over an order of magnitude speedup relative to known methods. Furthermore, the approach is parallel-structured so that it is suited for parallel implementation and further speedups. When used in conjunction with the recently developed local gravity approximations in conjunction with parallel computation, we anticipate MCPI will enable revolutionary speedups while ensuring accuracy.