I'm interested in numerically solving the geodesic equations of motion for a particle in the Schwarzschild geometry. Deriving the equations is no problem but I simply have no idea how to initialise and make sense of the output here. We have 4 spacetime coordinates given by $x^\alpha = (t(\lambda),r(\lambda),\theta(\lambda),\phi(\lambda))$, where $\lambda$ is an affine parameter and $t>0, 0<r<\infty, 0<\theta<\pi$ and $0<\varphi<2\pi$. Two of the equations of motion are given by \begin{equation} \frac{d \varphi}{d\lambda} = \frac{h}{r^2}, \\ \frac{dt}{d\lambda} = \frac{E}{\left(1-\frac{2GM}{c^2r}\right)}, \end{equation} where $E,h$ are the energy and the magnitude of the angular momentum. Has anyone references to numerically solving these equations? There is a Mathematica code which is very very cool but I'm more interested to do it myself as a learning exercise rather than just use this not to mention my Mathematica skills are almost as bad as my numerical.
For example, consider the second integral. After a certain number of integration steps the parameter $\lambda$ and $t$ should be different indicating time dilation. Again, I'm not really concerned with the physics of the problem more about the mathematical set-up and mainly some resources.
There are many different books on numerical solutions of ODEs and nearly all will consider coupled systems of ODEs, simply because any second- (or higher-) order ODE can be written as a system of first-order ODEs. A sample of some books I have used are:
As for methods, you could use finite differences, finite volumes, finite element, spectral methods - I don't believe you will be limited just because it is a system. What is more likely to be a problem is the nonlinearity in the geodesic equation!