A method is presented for solving boundary and initial value problems in celestial mechanics. In particular we consider the well-known Lambert TPBVP. The approach is quite general, however certain details in the transformed space boundary conditions pose challenges. We have been able to resolve these difficulties fully for the planar classical two-body problem, and we are engaged in a study to extend our numerical algorithm to the generally perturbed case. This method fuses three sets of ideas: (i) Picard Iteration, (ii) Orthogonal approximation, and notably, regularizing transformation of the equations of motion. Curiously, we find that a local-linearization-based shooting is not required, and we also illustrate that the method is not highly sensitive to the starting approximation. Two variants of the approach are considered, with the first model utilizing a Picard Iteration operating on the general differential equations in rectangular coordinates, which are approximated by Chebyshev polynomials. The second variant makes use of the KS transformation to render the unperturbed motion rigorously linear. These techniques combined improve the time interval over which the Picard Iteration converges, and increases the speed of convergence over all time intervals. A numerical study demonstrates excellent execution time efficiency, and shows that these algorithms are also attractive for parallelization if needed for further computational speedup. These new algorithms address improvements in the solutions of a fundamental problem in astrodynamics and should find widespread use in contemporary and future applications.