This paper is concerned with the numerical solution of the Euler-Lagrange equations of Mechanics, i.e. second-order ordinary differential equations where the variables satisfy some holonomic constraints. Typically such a set of equations in autonomous form can be written as \(M(x)x'' + F'(x)^ T z = G(x,x')\), together with \(F(x) = 0\), where \(F\) defines the constraints, \(M\) is the mass matrix (positive definite), the term \(G(x,x')\) is due to the external forces and \(z\) the Lagrange multipliers introduced by the constraint forces. The above system of equations, written as a first-order system in the variables \((x,x',z)\), is a differential-algebraic system of equations (DAE) of index three and standard techniques to solve it attempt to low the index and apply general numerical methods existing for DAEs of index 1 or 2. Here the method proposed by the authors follows the approach of analytical dynamics which consists in taking a suitable set of independent variables (Lagrangian coordinates) so that the constraints are automatically satisfied. Thus, in Sections 2 and 3 several smoothness assumptions are introduced in order to assure a local parametrization of the constraints manifold and to reduce the second-order equation to the manifold. In Section 4 some remarks about the numerical solution of the resulting second-order equation on the manifold are given. Finally, the numerical results with three test examples are presented and compared with those obtained with other solvers.