I'm currently trying to solve the geodesic equation in the most general setting possible. A curve $\gamma(s) = (\gamma^0(s), \gamma^1(s), \gamma^2(s), \gamma^3(s))$ satisfies the geodesic equation if
$$\frac{d^2\gamma^i}{ds^2} + \sum_{a,b}\Gamma_{ab}^i(\gamma(s))\frac{d\gamma^a}{ds}\frac{d\gamma^b}{ds} = 0$$ for $i\in\{0,1,2,3\}$. Making the substitution $\tau^i = d(\gamma^i)/ds$ you can reduce the system to one twice as big of the form $$\frac{d\tau^i}{ds} + \sum_{a,b}\Gamma_{ab}^i(\gamma(s))\frac{d\gamma^a}{ds}\frac{d\gamma^b}{ds} = 0$$ I tried to solve it as one usually does (that is, writing it in the form $y' = f(s, y)$) but, although factoring $\tau'$ is easy, one cannot do such a thing for $d\gamma^i/ds$: assuming one has the Levi-Civita connection (that is, $\Gamma^i_{ab} = \Gamma^i_{ba}$) one ends up with something of the form $$a(s) + b(s)\frac{d\gamma^i}{ds} + \Gamma_{ii}^i(\gamma(s))\left(\frac{d\gamma^i}{ds}\right)^2 = 0$$ Is there a method for solving one of these ODEs?; in summary, can one numerically solve an ODE of the form $(y')^2 + p(t)y' + q(t) = 0$?