Abstract
A sixth-order accurate scheme is presented for the solution of ODE systems supplemented by two-point boundary conditions. The proposed integration scheme is a linear multi-point method of sixth-order accuracy successfully used in fluid dynamics and implemented for the first time in astrodynamics applications. A discretization molecule made up of just four grid points attains a O(h 6) accuracy which is beyond the first Dahlquist's stability barrier. 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.