Solution (source code)

= Solution

Resolving the conserved energy into radial and azimuthal parts gives
$$
C=\frac12\left(\dot r^2+\frac{h^2}{r^2}\right)-\frac{\mu(1-\beta)}r.
$$
Using the release values of $C$ and $h$, and choosing the outward root,
$$
\boxed{v_r(r)=\dot r
=\sqrt{\frac{\mu}{a_p}\left[(2\beta-1)
-2(\beta-1)\frac{a_p}{r}-\left(\frac{a_p}{r}\right)^2\right]}}.
$$
Equivalently, with $x=r/a_p$,
$$
v_r=\sqrt{\frac{\mu}{a_p}}
\frac{\sqrt{(x-1)[(2\beta-1)x+1]}}x.
$$