How would I start to numerically compute trajectories of Kerr geodesics with constants of motion like in this wikipedia page. I want to recreate trajectories like in this picture in Matlab.
A seminal paper on Kerr Geodesics is Wilkins. The necessary equations are found at 2 and 3.
Note that it is not trivial to implement these equations by plugging into an RK4 integrator because of the square roots in the R and Theta "potential" functions.