I am looking for an algorithm to solve numerically the following system of differential equations. It is the geodesic equation on for example a 2-dimensional manifold. So, let $\gamma(t)=(\gamma_1(t),\gamma_2(t))$ be the geodesic passing through point $\gamma(0)=(x_0,y_0)$ with initial velocity $\gamma'(0)=(v_1,v_2)$. I know Runge-Kutta method for solving a differential equation of first order but no idea about a system of second order differential equations.
Eqaution of geodesic: $$ \gamma''_1(t)+ \Gamma_{111}(\gamma(t)) \hspace{0.2cm}\gamma_1'(t)\gamma'_1(t)+\Gamma_{121}(\gamma(t)) \hspace{0.2cm}\gamma_2'(t)\gamma'_1(t)+\Gamma_{112}(\gamma(t)) \hspace{0.2cm}\gamma_1'(t)\gamma'_2(t)+ \Gamma_{122}(\gamma(t)) \hspace{0.2cm}\gamma_2'(t)\gamma'_2(t)=0\\ \gamma''_2(t)+ \Gamma_{211}(\gamma(t)) \hspace{0.2cm}\gamma_1'(t)\gamma'_1(t)+\Gamma_{221}(\gamma(t)) \hspace{0.2cm}\gamma_2'(t)\gamma'_1(t)+\Gamma_{212}(\gamma(t)) \hspace{0.2cm}\gamma_1'(t)\gamma'_2(t)+ \Gamma_{222}(\gamma(t)) \hspace{0.2cm}\gamma_2'(t)\gamma'_2(t)=0$$
The set of functions $\Gamma_{ijk}$ are Christoffel symbols and are real valued functions of 2-variables.
I got a suggestion about reducing the order of equations by defining extra functions $f=\gamma'_1$ and $g=\gamma'_2$ and writing down the associated relations and forming a new system. Yet I have my question unanswered because I don't know the algorithm even I found a similar question here. But the mentioned question tries to give an algorithm to find two functions $z,y$.
Edited question:
Can any one generalize mentioned question's algorithm for my problem?