Computing Trajectory Equations of Kerr Geodesics

I want to numerically solve the trajectory equations of a Kerr geodesic given by wikipedia in Matlab.

The trajectories look like: I implemented the equations and solved it with the standard Runge-Kutta method. Here it is my code:

% Initial parameter (r, theta, phi, t)
r       = x(1);
theta   = x(2);
phi     = x(3);
t       = x(4);

% Constants of motion
mu = const(1);
E = const(2);
L = const(3);
Q = const(4);

% Abbreviations
Sigma = r^2 + a^2 * cos(theta)^2;
Delta = r^2 - 2 * M * r + a^2;
K = Q + (L - a*E)^2;
X = a^2*(E^2 + mu) - L^2 - Q;
R = (E^2 + mu)*r^4 - 2*M*mu*r^3 + X*r^2 + 2*M*K*r -a^2*Q;
O = Q + cos(theta)^2*(a^2*(E^2 + mu) - L^2/sin(theta)^2);

% ODE from O Neill Chapter 4.2.2 Geodesics
dx(1) = sqrt(R / Sigma^2);
dx(2) = sqrt(O / Sigma^2);
dx(3) = ((L - a*sin(theta)^2*E) / sin(theta)^2 + a*((r^2 + a^2)*E - L*a)/(Delta) ) / Sigma^2;
dx(4) = (a*(L - a*E*sin(theta)^2) + (r^2 + a^2)*((r^2 + a^2)*E - L*a)/(Delta)) / Sigma^2;

I am a little confused by the R- and Q-term because there is a square root. How can I handle the signs of these in the calculations?

I know there exists an alternative form of the equations with the components of momenta, but I do not want to use this approach.

• What is the meaning of the sign on the equations? – nicoguaro May 31 at 16:58
• @nicoguaro By the signs I mean the "+" and "-" of the R- and Q-term, because there is the squareroot. – almost Jun 2 at 6:47
• I know. But that does not answer my question. To know how to handle them numerically it is good idea to know what their role is. – nicoguaro Jun 2 at 13:52
• @nicoguaro as per the original question on Physics SE, they would be the "potential" functions. Also, in the answer there, a very useful Wilkins paper has been referenced. – Anton Menshov Jun 5 at 5:18
• @nicoguaro don't know that yet. Plan to read the paper some time later, as this question is a bit relevant (very far-fetched) to a side project of mine. – Anton Menshov Jun 6 at 19:41