  In this work, we explore the eternal collapsing phenomenon of a stellar
system (e.g., a star) within the framework of $f(R)$ gravity and investigate
some new aspects of the continued homogeneous gravitational collapse with
perfect fluid distribution. The exact solutions of field equations have been
obtained in an independent way by the parameterization of the expansion scalar
($\Theta$) governed by the interior spherically symmetric FLRW metric. We
impose the Darmois junction condition required for the smooth matching of the
interior region to the Schwarzschild exterior metric across the boundary
hypersurface of the star. The junction conditions demand that the pressure is
non-vanishing at the boundary and is proportional to the non-linear terms of
$f(R)$ gravity, and the mass function $m(t, r)$ is equal to Schwarzschild mass
$M$. The eight massive stars, namely $Westerhout 49-2, BAT99-98, R136a1,
R136a2, WR 24, Pismis 24-1$, $\lambda- Cephei$, and $\beta -Canis Majoris$ with
their known astrophysical data (masses and radii) are used to estimate the
numerical values of the model parameters which allows us to study the solutions
numerically and graphically. Here we have discussed two $f(R)$ gravity models
describing the collapse phenomenon. The singularity analysis of models is
discussed via the apparent horizon and we have shown that stars tend to
collapse for an infinite co-moving time in order to attain the singularity (an
eternal collapsing phenomenon). We have also shown that our models satisfy the
energy conditions and stability requirements for stellar systems.
