I have these 5 equations: \begin{align} \dot\psi&=\frac{p_\psi}{2r} -\frac{\sin(\psi)}{r} \\ \dot r&=\cos(\psi) \\ \dot h&=\sin(\psi) \\ \dot p_\psi&=\left(\frac{p_\psi}{r} -p_h\right)\cos(\psi)+\left(\frac{2\sigma r}{a^2} +p_r\right)\sin(\psi) \\ \dot p_r&=\frac{p_\psi}{r}\left(\frac{p_\psi}{4r} -\frac{\sin(\psi)}{r}\right)+\frac{2\sigma}{a^2}\big(1-\cos(\psi)\big) \\ \end{align}
Where $ \psi , h, r, p_\psi, p_r$ are functions of $s$. with five boundary conditions at $s=0$ $$r(0),h(0),\psi(0),p_\psi(0),p_r(0)$$
I already can solve them using software like Mathematica and Maple, but I want to write a code to solve them numerically, and I don't know which method can be used. Do you have any suggestions?