Solution (source code)

= Solution

Spherical symmetry puts any geodesic in a plane through the origin, chosen as $\theta=\pi/2$. The cyclic coordinates give
$$
E=(1-r_s/r)\dot t,\qquad L=r^2\dot\phi.
$$
The null condition gives
$$
\dot r^2=F(r)=E^2-(1-r_s/r)\frac{L^2}{r^2}.
$$
For $u=1/r$, $(du/d\phi)^2=\dot r^2/L^2$, hence
$$
\boxed{\left(\frac{du}{d\phi}\right)^2=\frac1{b^2}-u^2+r_su^3,\qquad b=L/E.}
$$