Runge-Kutta Approach to Relativistic Kepler Problem

I am attempting to solve the relativistic Kepler problem using fourth order Runge-Kutta, however I am having trouble and am unsure why I simply cannot get the correct orbits no matter how I adjust the initial values. My orbits are not even incorrect but they are nonsensical, I was wondering if there are issues with the Runge-Kutta algorithm when solving such problems? The geodesic equations in Schwarzschild spacetime (in geometric units) are:

begin{align}
&frac{d^2r}{dtau^2} = frac{M}{r(r – 2M)} left(frac{dr}{dtau}right)^2-frac{M(r – 2M)}{r^3} left(frac{dt}{dtau}right)^2 + (r – 2M) left(frac{dphi}{dtau}right)^2
[0.5em]
&frac{d^2t}{dtau^2} = – frac{2M}{r(r – 2M)} frac{dr}{dtau} frac{dt}{dtau}
[0.5em]
&frac{d^2phi}{dtau^2} = – frac{2}{r} frac{dr}{dtau} frac{dphi}{dtau}
end{align}

We can convert the above to a system of 6 1st order ODE’s by setting:
begin{align}
y_1 = r [0.5em]
y_2 = frac{dr}{dtau} [0.5em]
y_3 = phi [0.5em]
y_4 = frac{dphi}{dtau} [0.5em]
y_5 = t [0.5em]
y_6 = frac{dt}{dtau}
end{align}

So then, we obtain 6 1st order ODEs which can be solved using the fourth order Runge – Kutta. For example, I tested my code on mercury with the following parameters (geometrized units):
begin{align}
&text{Mass of sun: } 1474 m [0.5em]
&text{Mecury radius at perihelion: } 46.002times 10^9 m
end{align}

I am not sure about the initial values for the Runge Kutta so I tried a range such as
$$
[y_1(0), y_2(0), y_3(0), y_4(0), y_5(0), y_6(0)] = [46.002times 10^9, 0, 0, 1times 10^{-16}, 0, 1]
$$
Despite the range of initial values i have tried, I cannot seem to get even a decent orbit.

Is there a problem with using Runge-Kutta to solve these equations or should there be a more deterministic way to obtain the initial values?

Mmercury = 3.3010 10^23; RAmercury = 5.7909227 10^10; RPmercury = 4.6002 10^10; Emercury = 0.20563593; Pmercury = 87.97; vp=58983;

c0 = 299792458.; GM = 1.32712 10^(20); mS = GM/c0^2; M=mS/RPmercury ;