This is a follow-up question to the answer given at What is the exact gravitational force between two masses including relativistic effects?. Unfortunately the author hasn't been online for a few years and therefore does not answer comments any more.
In the answer given there the differential equation of motion in Schwarzschild coordinates was
$$\ddot{r} = -\frac{G m}{r^2} + r\dot{\theta}^2 - \frac{ {\color{red} 3} G m}{c^2}\dot{\theta}^2$$
for the radial acceleration and
$$\ddot{\theta} = -\frac{2}{r}\dot{r}\dot{\theta}$$
for the angular acceleration. When I plot the path for an object near the speed of light, with this formula I get a stable orbit at $r_0=2 r_s$:
But shouldn't that be at $r_0=1.5 r_s$, the photon sphere? With that formula the orbiting particle would fall into the black hole very quickly, for example, with $v_0=0.999c$ at $r_0=1.6 r_s$:
When I replace the term 3Gm/c² with 2GM/c² so that
$$\ddot{r} = -\frac{G m}{r^2} + r\dot{\theta}^2 - \frac{ {\color{red} 2} G m}{c^2}\dot{\theta}^2$$
I get the expected result with a stable orbit right at the photon sphere (initial velocity again $v_0=0.999c$):
So my question is: is the formula wrong and the factor 3 needs to be replaced with a factor of 2, or are there different minimum-distances for stable orbits, one for particles, and one for photons? Or did I miss something else? Wikipedia says:
The radius of the photon sphere, which is also the lower bound for any stable orbit, is $1.5 r_s$
so I would expect that also particles with mass should stay in orbit if they are close to the speed of light and slightly above the photon sphere.
For reproduction of the problem the Mathematica-code as I believe it to be correct is available (with the factor 2 instead of 3)