  We use Lagrangian formalism and jet spaces to derive a computational model to
simulate multibody dynamics with holonomic constraints. Our approach avoids the
traditional problems of drift-off and spurious oscillations. Hence even long
simulations remain physically relevant. We illustrate our method with several
numerical examples.
