Is there something wrong with this numerical simulation of Schwarzschild photon orbits?

There seem to be several confusions here. Massive and massless particles behave qualitatively differently, even if the massive particle is traveling very fast.

• The minimum radius for a stable orbit for a massive particle is $3 r_s$. Circular orbits above this radius are all stable.
• Massless particles only have circular orbits at the photon sphere, $(3/2) r_s$. These orbits are not stable, Wikipedia is wrong. Massless particles also obey a different equation of motion.

The other confusion is that what your simulations are showing has nothing to do with stability. Your particles are falling into the center because you're not giving them the right initial velocity. It's just like classical mechanics: if you suddenly took away half the Earth's speed, it would start falling inward. In order to initialize them to the right initial velocity, you need to solve for $\dot{\theta}$ so that $\ddot{r} = 0$.

This is in contrast to the massless case, where the initial velocity is already determined for you (i.e. it's the speed of light).

Thanks to the hint given by knzhou I figured out that if one wants to give the particle a proper initial velocity of $v_0$, the initial velocity in terms of Schwarzschild coordinates $v_i$ would then be

$$\dot{\theta}(0)\cdot r(0) = \frac{{v\perp}_0}{ \color{green}{\sqrt{ 1-v_0^2/c^2}}}$$

for the transversal component, and

$$\dot{r}(0) = \frac{{v\parallel}_0\cdot \color{blue}{\sqrt{1-r_s/r_0}}}{ \color{green}{\sqrt{ 1-v_0^2/c^2}}}$$

for the radial component since one has to compensate for the gravitational length contraction (blue) and the length contraction due to the particle's velocity (green) with respect to the particle's proper time.

Now I get the expected results: a circular orbit with transversal initial velocity around the photon sphere, and a stationary particle with outwards initial velocity around the event horizon when v0 is set close to c.