Abstract
The present paper presents an accurate scheme for the solution of boundary value problems with two-point nonlinear boundary conditions. The proposed scheme is a linear multi-point method of sixth-order accuracy successfully used in uid dynamics and here implemented for the rst time in astrodynamics applications. It is an optimal scheme since a discretization molecule made up of just four grid points assures an h6 order of accuracy. This kind of discretization allows to attain an accuracy beyond the rst Dahlquist's stability barrier and simultaneously has a simple formulation and numerical e ciency. Astrodynamics applications concern the computation of libration point halo orbits, in the restricted three- and four-body models, and the design of an optimal control strategy for a low thrust libration point mission.