- Modified Chebyshev-Picard iteration methods are presented for solving boundary value problems. Chebyshev polynomials are used to approximate the state trajectory in Picard iterations, while the boundary conditions are maintained by constraining the coefficients of the Chebyshev polynomials. Using Picard iteration and Clenshaw-Curtis quadrature, the .presented methods iteratively refine an orthogonal function approximation of the entire state trajectory, in contrast to step-wise, forward integration approaches, which render the methods well-suited for parallel computation because computation of force functions along each path iteration can be rigorously distributed over many parallel cores with negligible cross communication needed. The presented methods solve optimal control problems through Pontryagin's principle without requiring shooting methods or gradient information. The methods are demonstrated to be computationally efficient and strikingly accurate when compared with Battin's method for a classical Lambert's problem and with a Chebyshev pseudospectral method for an optimal trajectory design problem. The reported simulation results obtained on a serial machine suggest a strong basis for optimism of using the presented methods for solving more challenging boundary value problems, especially when highly parallel architectures are fully exploited.