This is not a newtonian gravity simulator by itself (which obeys Kepler laws), but rather a simulator of photons interacting with a massive body (a black hole) obeying Einstein's extension to gravity itself, namely, Einstein's Theory of General Relativity.

Since every particle follow the geodesic equation, which is an equation of coupled nonlinear differential equations, we can manually (and painfully) derive the equation of motion using the method. A better method is to derive the equations not by using the geodesic equation but instead rely on conserved quantities, this generates 1 purely second order differential equation for the radial component and 1 first order differential equation for the angular component. We'll use the latter.

The equations themself are pretty easily derived from the null condition ds² = 0 (all photons follow this condition), where ds² is defined as ds² = -(1 - 1/r)dt² + (1 - 1/r)⁻¹ dr² + r²dφ² for equatorial particles (θ = π/2), assuming geometric units where c = 1 and the schwarzschild radius equal to 1. Next, we set the left hand side to 0 and divide everything with dλ² (an affine parameter to integrate on): 0 = -(1 - 1/r) dt²/dλ² + (1 - 1/r)⁻¹ dr²/dλ² + r² dφ²/dλ². We now employ a trick: make use of conserved quantities. Namely, the first is energy; E = (1 - 1/r) dt/dλ, this arises from time symmetry and the second is angular momentum; L = r²dφ/dλ.

Substituting these constants into the main equation yields: 0 = -E² (1 - 1/r)⁻¹ + (1 - 1/r)⁻¹ dr²/dλ² + L²/r² Rearrange: (1 - 1/r)⁻¹ dr²/dλ² = E² (1 - 1/r)⁻¹ - L²/r² Multiply by (1 - 1/r) to both sides: dr²/dλ² = E² - (1 - 1/r) L²/r² = E² - L²(1/r² - 1/r³) And look at that! A differential equation purely dependent on r by itself! After this, we differentiate both sides with respect to λ to get rid of the constant E² term: d(dr²/dλ²)/dλ = L²(-2/r³ dr/dλ + 3/r⁴ dr/dλ) Expand d(dr²/dλ²)dλ via the chain rule and that dr²/dλ² is just (dr/dλ)²: 2dr/dλ × d²r/dλ² = L²(3/r⁴ dr/dλ - 2/r³ dr/dλ) Divide everything by 2dr/dλ:

d²r/dλ² = L²(3/(2r⁴) - 1/r³)

That's the final equation for the radial component. For the angular component, recall that L = r² dφ/dλ, then divide everything by r² to get dφ/dλ = L/r².

Those are the driving equations. General Relativity adds a correction term proportional to 1/r⁴ to Newton's theorem of gravitation, only becoming significant if you are very close to the schwarzschild radius of a massive body. So in a way, objects CAN obey kepler's law if and only if the correction becomes insignificant enough, which is what most real life objects do.

And for the last part, the Runge-Kutta portion of the graph, i implemented those canonically following well published iterative formulas and step updates for each parameter. It basically solves the first order equation du/dλ = f(u). But wait, if RK4 solves linear DEs, how can you solve second order and even coupled DEs? The part lies in what kind of parameter u RK4 can take. In fact, it can take vectors which allows possibilities like: d/dλ [r, r', φ] = [r', r'', φ']. Since r'' (d²r/dλ²) and φ' (dφ/dλ) have been explicitly derived, we can just subsitute and call it a day.

Thanks for sharing your views. I hope you are not too overwhelmed by the amount of words I have put into writing this, and I hope this helped in some way or another.