Setting up the Hamiltonian of the system, then substituting into the equations and evaluating the partial derivatives and simplifying, a set of coupled 2 N first order ODEs in the coordinates q i and momenta p i are obtained.