Solution (source code)

= Solution

The radial first integral for an equatorial <null geodesic> in <Schwarzschild spacetime> is
$$
\dot r^2+\left(1-\frac{2M}{r}\right)\frac{L^2}{r^2}=E^2.
$$
Since $L=r^2\dot\phi$, regarding $r$ as a function of $\phi$ gives
$$
\dot r=\frac{dr}{d\phi}\dot\phi
=\frac{L}{r^2}\frac{dr}{d\phi}.
$$
Set $y=1/r$. Then $dr/d\phi=-y'/y^2$ and therefore $\dot r=-Ly'$. Division by $L^2$ gives the <equatorial null orbit in Schwarzschild spacetime>
$$
\boxed{(y')^2+y^2-2My^3=\frac{E^2}{L^2}=\frac1{b^2},}
$$
where $b=L/E$ is the <impact parameter>. Differentiating with respect to $\phi$ gives
$$
2y'(y''+y-3My^2)=0.
$$
By continuity this yields
$$
\boxed{y''+y=3My^2.}
$$