The solution of initial value problems (IVPs) provides the evolution of dynamic system state history for given initial conditions. Solving boundary value problems (BVPs) requires finding the system behavior where elements of the states are defined at different times. This dissertation presents a unified framework that applies modified Chebyshev-Picard iteration (MCPI) methods for solving both IVPs and BVPs. Existing methods for solving IVPs and BVPs have not been very successful in exploiting parallel computation architectures. One important reason is that 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. The proposed MCPI methods are inherently parallel algorithms. Using Chebyshev polynomials, it is straightforward to distribute the computation of force functions and polynomial coefficients to different processors. Combining Chebyshev polynomials with Picard iteration, MCPI methods iteratively refine estimates of the solutions until the iteration converges. The developed vector-matrix form makes MCPI methods computationally efficient. The power of MCPI methods for solving IVPs is illustrated through a small perturbation from the sinusoid motion problem and satellite motion propagation problems. Compared with a Runge-Kutta 4-5 forward integration method implemented in MATLAB, MCPI methods generate solutions with better accuracy as well as orders of magnitude speedups, prior to parallel implementation. Modifying the algorithm to do double integration for second order systems, and using orthogonal polynomials to approximate position states lead to additional speedups. Finally, introducing perturbation motions relative to a reference motion results in further speedups. The advantages of using MCPI methods to solve BVPs are demonstrated by addressing the classical Lambert’s problem and an optimal trajectory design problem. MCPI methods generate solutions that satisfy both dynamic equation constraints and boundary conditions with high accuracy. Although the convergence of MCPI methods in solving BVPs is not guaranteed, using the proposed nonlinear transformations, linearization approach, or correction control methods enlarge the convergence domain. Parallel realization of MCPI methods is implemented using a graphics card that provides a parallel computation architecture. The benefit from the parallel implementation is demonstrated using several example problems. Larger speedups are achieved when either force functions become more complicated or higher order polynomials are used to approximate the solutions.