We present an efficient variational integrator for simulating multibody systems. Variational integrators reformulate the equations of motion for multibody systems as discrete Euler-Lagrange (DEL) equation, transforming forward integration into a root-finding problem for the DEL equation. Variational integrators have been shown to be more robust and accurate in preserving fundamental properties of systems, such as momentum and energy, than many frequently used numerical integrators. However, state-of-the-art algorithms suffer from O(n 3 ) complexity, which is prohibitive for articulated multibody systems with a large number of degrees of freedom, n, in generalized coordinates. Our key contribution is to derive a quasi-Newton algorithm that solves the root-finding problem for the DEL equation in O(n), which scales up well for complex multibody systems such as humanoid robots. Our key insight is that the evaluation of DEL equation can be cast into a discrete inverse dynamic problem while the approximation of inverse Jacobian can be cast into a continuous forward dynamic problem. Inspired by Recursive Newton-Euler Algorithm (RNEA) and Articulated Body Algorithm (ABA), we formulate the DEL equation individually for each body rather than for the entire system, such that both inverse and forward dynamic problems can be solved efficiently in O(n). We demonstrate scalability and efficiency of the variational integrator through several case studies.