How to solve equation of motion for Orbit in Einsteinian Mechanics? Hi i am trying to solve the following equation using runge-kutta method but it doesn't seem to work. $$\frac{dr}{d\phi}=\sqrt{\frac{2Er^4}{l^2} + \frac{2GM}{l^2} -r^2 + 2GMr}.$$  I have written a python script to solve this using runge-kutta method. My question is : is there any other way to solve this equation? Also will this have an analytical solution. I want $r$ in terms of $\phi$.
PS:
$E,M,l,G$ are constants.
 A: The problem is probably with the turning points $dr/d\phi = 0$ (also notice that you have to switch the sign of the right hand side after passing through it!), since if you overshoot them even by a little, you will get an imaginary value for the derivative. A solution for this problem is one of the following:

*

*To work with a second-order equation for $d^2r/d\phi^2 = f(r,E,l,dr/d\phi)$, where all square-root expressions are replaced by $dr/d\phi$.

*To choose a parametrization $r(\xi) = (r_1 + r_2) + (r_1-r_2)\sin \xi$, where $r_1,r_2$ are the roots of the equation $dr/d\phi = 0$. You should then be able to re-express the equation as $d \xi/d\phi = g(\xi,E,l)$ without the turning-point singularities.

*Work with the implicit equation
$$\int_{r_0}^{r} \frac{dr'}{\pm\sqrt{2Er'^2/l^2 + ...}} = \phi-\phi_0$$
You can integrate the left-hand side by any method for definite integrals, and then you can get the series of values $r,\phi$ along the trajectory (but remember, you have to switch the sign/reflect the direction of integration once encountering the turning point).

*You could just use your original algorithm and put checks on whether you get an imaginary value for $dr/d\phi$. If you do, you can roll back a step (which you have to have saved), and either make the step smaller, or if $dr/d\phi$ is small enough, switch the sign of $dr/d\phi$ to start integrating away from the turning point (depending how you do this, this may be associated with increased numerical error).

