  We study the Starobinsky or $R^2$ model of $f(R)=R+\alpha R^2$ for neutron
stars with the structure equations represented by the coupled differential
equations and the \emph{polytropic} type of the matter equation of state. The
junction conditions of $f(R)$ gravity are used as the boundary conditions to
match the Schwarschild solution at the surface of the star. Based on these the
conditions, we demonstrate that the coupled differential equations can be
solved \emph{directly}. In particular, from the dimensionless equation of state
$\bar{\rho} = \bar{k}\, \bar{p}^{\,\gamma}$ with $\bar{k}\sim5.0$ and
$\gamma\sim0.75$ and the constraint of $\alpha\lesssim {1.47722}\times 10^{7}\,
\text{m}^2$, we obtain the \emph{minimal} mass of the NS to be around 1.44
$M_{\odot}$. In addition, if $\bar{k}$ is larger than 5.0, the mass and radius
of the NS would be smaller.
