I have a particular Riemannian metric $g$ defined on a subset of Euclidean space $E \subset \mathbb{R}^n$. Although this is beyond the question, I have $E = \{x:\|x\|_2 > 1\}$. So really, $g$ is a smooth map from $E$ to the space of positive definite matrices. Excitingly, I proved the resulting Riemannian manifold $(g,E)$ is geodesically complete. As a result, I want to numerically compute the geodesics connecting pairs of points in $E$.
As you know, a geodesic on my Riemannian manifold $(E,g)$ is a solution $x(t)$ to the following ODE: $$\ddot{x}^k(t)+\Gamma_{ij}^k(x(t))\dot{x}^i(t) \dot{x}^j(t)=0,$$ where $\Gamma_{ij}^k := \frac{1}{2}g^{il}(g_{kl,j} + g_{jl,k}-g_{jk,l})$. Here $g^{ij}$ is the $ij$th entry of the inverse of $g$ and $g_{ij,k}:=\frac{\partial g_{ij}}{\partial x^k}$.
The smooth manifold structure is simply the submanifold structure induced from the ambient Euclidean space.
I tried implementing the geodesic equations in Matlab and solving them in the usual way (ODE45). However, I found that the ODE is quite unstable. For anything other than small values of $t$, the numerical solution goes unstable. I think the ODE is stiff. I found another solve that works somewhat better, but only for solving $x(t)$ for small values $0 \leq t \leq T$.
I have the following 3 questions:
What is the best integrator for solving the geodesic equation? Excuse my limited knowledge of the topic, but it seems variational or symplectic integrators might be of use. What is the best integrator? Clearly it is not ODE45 or the Euler method.
What is the best way to solve geodesic equation for 2-point boundary equations?
Are there any pre-written libraries dedicated to this stuff? So far, I cannot find any. I have come across some for General relativity computations, but clearly that is not the case for me.