I am trying to solve the following set of differential equations for plotting particle trajectories (using 4th order Runge-Kutta method): $$\frac{dr}{dt}=-\Delta\left[\frac{2}{r}\left(1+\frac{a^2}{r^2}\right)\right]^{1/2}\left(r^2+a^2+\frac{2a^2}{r}\right)^{-1}$$ $$\frac{d\phi}{dt}=\frac{2a^2}{r}\left(r^2+a^2+\frac{2a^2}{r}\right)^{-1}$$ where $\Delta=r^2+a^2-2r$ and $a$ is a constant. We require to plot the trajectory in the $r-\phi$ plane (or equivalently, in the $x-y$ plane with $x=\sqrt{r^2+a^2}\cos\phi$ and $y=\sqrt{r^2+a^2}\sin\phi$).
The above equations can also be written as a single differential equation as follows: $$\frac{d\phi}{dr}=-\frac{2a}{r\Delta}\left[\frac{2}{r}\left(1+\frac{a^2}{r^2}\right)\right]^{-1/2}$$
The difficulty I am facing in solving the problem is that when $\Delta=0$ (called horizons), $\dfrac{dr}{dt}$ and $\dfrac{dr}{d\phi}$ becomes zero, but $\dfrac{d\phi}{dt}$ remains finite there. For this reason, the trajectory has discontinuities at two values of $r$ when $\Delta=0$. The following is a plot of the trajectories for $a=0.8$ and the initial conditions $(r,\phi)=(8,0)$ at $t=0$:
In the figure, the axes are $x=\sqrt{r^2+a^2}\cos\phi$, $y=\sqrt{r^2+a^2}\sin\phi$, $\mu=1$ (already used in the equations) and the dashed lines correspond to the horizons ($\Delta=0$). Note that the behavior of the trajectory changes at the horizon.
It would be very helpful if someone could suggest me how to deal with the discontinuities (shown in dashed lines) in the Runge-Kutta algorithm for the particle trajectories.
