Solution (source code)

= Solution

The radial equation is
$$
\ddot r=-\frac{GM}{r^2}+\frac{\Lambda c^2}{3}r.
$$
Multiplying by $\dot r$ and integrating gives
$$
\frac12\dot r^2-\frac{GM}{r}
-\frac{\Lambda c^2}{6}r^2=E.
$$
Since $\rho$ is an energy density,
$$
M=\frac{4\pi r^3}{3c^2}\rho.
$$
Set $r=ax_0$, divide the energy equation by $r^2/2$, and define
$$
\frac{2E}{x_0^2}=-kc^2.
$$
Then
$$
\boxed{
\left(\frac{\dot a}{a}\right)^2
=\frac{8\pi G}{3c^2}\rho-\frac{kc^2}{a^2}
+\frac{\Lambda c^2}{3}}.
$$