I have encountered the following system of differential equations in lagrangian mechanics. Can you suggest a numerical method, with relevant links and references on how can I solve it, and the implementation in C (if possible) Also, is there a shorter implementation on Matlab or Mathematica?
\begin{align*} mx \dot y^2 + mg\cos(y) - Mg - (m+M)\,\ddot x &= 0 \\ g\sin(y) + 2\dot x\dot y + x \,\ddot y &= 0 \end{align*}
where $\dot x$ or $\dot y$ are time derivatives, and the double dots indicate a 2nd derivative wrt time.