Solution (source code)

= Solution

For an infinitesimal test mass $\delta m$ on the surface of a homogeneous sphere, its Newtonian <energy>[mechanical energy] is
$$
E=\frac12\delta m\,\dot R^2-\frac{GM\delta m}{R}.
$$
The spatially flat case corresponds to zero specific energy. Using
$$
M=\frac{4\pi}{3}\rho R^3
$$
and dividing by $\delta mR^2/2$ gives
$$
\left(\frac{\dot R}{R}\right)^2=\frac{8\pi G}{3}\rho.
$$
Since $R$ is proportional to the <scale factor> $a$, this is the <Newtonian derivation of the flat Friedmann equation>
$$
\boxed{H^2=\left(\frac{\dot a}{a}\right)^2=\frac{8\pi G}{3}\rho}.
$$

Now take $a(t_0)=1$ and $\rho=\rho_0a^{-4}$. Put
$$
C=\sqrt{\frac{8\pi G\rho_0}{3}}.
$$
The contracting branch of the <Friedmann equation> is
$$
\frac{\dot a}{a}=-Ca^{-2},
\qquad
\dot a=-\frac Ca.
$$
Integrating $a\,da=-C\,dt$ from $t_0$ gives
$$
\frac12(a(t)^2-1)=-C(t-t_0).
$$
Therefore
$$
\boxed{a(t)=\sqrt{1-2C(t-t_0)}}
$$
on the contracting interval, reaching $a=0$ at $t=t_0+1/(2C)$.