The inverse dynamics of flexible multibody systems is formulated as a two-point boundary value problem for an index-3 differential-algebraic equation (DAE). This DAE represents the equation of motion with kinematic and trajectory constraints. For so-called nonminimum phase systems, the remaining dynamics of the inverse model is unstable. Therefore, boundary conditions are imposed not only at the initial time but also at the final time in order to obtain a bounded solution of the inverse model. The numerical solution strategy is based on a reformulation of the DAE in index-2 form and a multiple shooting algorithm, which is known for its robustness and its ability to solve unstable problems. The paper also describes the time integration and sensitivity analysis methods that are used in each shooting phase. The proposed approach does not require a reformulation of the problem in input-output normal form, which is known from nonlinear control theory. It can deal with serial and parallel kinematic topology, minimum phase and nonminimum phase systems, and rigid and flexible mechanisms.