We discuss both the theoretical framework and implementation issues for the simulation of holonomic mechanical systems using automatic differentiation. Our approach allows the direct simulation of such systems based on its Lagrangian. The Lagrangian may be given as explicit equations or as an algorithm containing control structures such as loops. The method is illustrated by the well-known ball and beam system.