We study the Starobinsky or $R^2$ model of $f(R)=R+α 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}^{\,γ}$ with $\bar{k}∼5.0$ and $γ∼0.75$ and the constraint of $α\lesssim {1.47722}\times 10^{7}\, \text{m}^2$, we obtain the \emph{minimal} mass of the NS to be around 1.44 $M_{⊙}$. In addition, if $\bar{k}$ is larger than 5.0, the mass and radius of the NS would be smaller. ; Comment: 17 pages, 3 figures