  Using numerical methods, we examine, under a Gowdy symmetry assumption, the
dynamics of nonlinearly perturbed FLRW fluid solutions of the
Einstein-Euler-scalar field equations in the contracting direction for linear
equations of state $p = K\rho$ and sound speeds $0\leq K<1/3$. This article
builds upon the numerical work from \cite{BMO:2023} in which perturbations of
FLRW solutions to the Einstein-Euler equations with positive cosmological
constant in the expanding time direction were studied. The numerical results
presented here confirm that the instabilities observed in
\cite{BMO:2023,MarshallOliynyk:2022} for $1/3<K<1$, first conjectured to occur
in the expanding direction by Rendall in \cite{Rendall:2004}, are also present
in the contracting direction over the complementary parameter range $0\leq
K<1/3$. Our numerical solutions show that the fractional density gradient of
the nonlinear perturbations develop steep gradients near a finite number of
spatial points and become unbounded towards the big bang. This behaviour, and
in particular the characteristic profile of the fractional density gradient
near the big bang, is strikingly similar to what was observed in the expanding
direction near timelike infinity in the article \cite{BMO:2023}.
