<p>The solution of two-point boundary value problems (TPBVPs) in multibody system dynamics is a tough challenge. In this article, a double shooting method (DSM) tailored for TPBVPs is presented. The method builds upon the classical single shooting method while incorporating optimization strategies to overcome numerical instability and sensitivity to initial guesses. Conventional shooting methods are associated with some numerical problems: The numerical gradient computation for updating initial values causes problems in many cases, since the disturbance parameter for determining the numerical derivative must be selected appropriately in order to achieve sufficient accuracy. Moreover, integrating the initial value problem can be numerically challenging, especially when the interval is large and the differential equations have unstable modes. The DSM is designed to address these challenges. Therefore, the method incorporates a discrete costate variable approach to cope with the numerical gradient computation, which is often the Achilles’ heel of classical shooting methods. Finally, the method is formulated using a discrete implicit integration scheme to determine, as a first example, a double pendulum upswing maneuver, as a second example, the energy-optimal control of a robot multibody model and, as a third example, trajectories of the restricted three-body problem.</p>