We present a new algorithm for approximating solutions of two-point boundary value problems and prove theorems that give conditions under which the solution must exist and the algorithm generate approximations that converge to it. We show how to make the algorithm computationally efficient and demonstrate how the full method works both when guaranteed to do so and more broadly. We demonstrate that the method compares well against other methods commonly used in the literature. We also prove a theorem on existence of solutions of certain multi-dimensional Volterra integral equations and use it to show that the Parker-Sochacki method of introducing auxiliary variables, used to make the new algorithm computationally efficient, can be effectively applied to these Volterra integral equations in order to approximate their solutions by means of a Picard iteration scheme. Finally, we extend the existence theorem for solutions of two-point boundary valueproblems and prove that the new algorithm can be modified to approximate solutions in this case.
