  We discuss the Mass -Radius diagram for static neutron star models obtained
by the numerical solution of modified Tolman-Oppenheimer-Volkoff equations in
$f(R)$ gravity where the Lagrangians $f(R)=R+\alpha R^2 (1+\gamma R)$ and
$f(R)=R^{1+\epsilon }$ are adopted. Unlike the case of the perturbative
approach previously reported, the solutions are constrained by the presence of
an extra degree of freedom, coming from the trace of the field equations. In
particular, the stiffness of the equation of state determines an upper limit on
the central density $\rho_c$ above which the the positivity condition of
energy-matter tensor trace $T^{\rm m}=\rho - 3 p$ holds. In the case of
quadratic f(R)-gravity, we find higher masses and radii at lower central
densities with an inversion of the behavior around a pivoting $\rho_c$ which
depends on the choice of the equation of state. When considering the cubic
corrections, we find solutions converging to the required asymptotic behavior
of flat metric only for $\gamma < 0$. A similar analysis is performed for
$f(R)=R^{1+\epsilon }$ considering $\epsilon$ as the leading parameter. We work
strictly in the Jordan frame in order to consider matter minimally coupled with
respect to geometry. This fact allows us to avoid ambiguities that could emerge
in adopting the Einstein frame.
