  In the present paper we investigate non-perturbatively and self-consistently
the structure of neutron stars in $R$-squared gravity by simultaneously solving
the interior and exterior problem. The mass-radius relations are obtained for
several equations of state and for wide range of the $R$-squared gravity
parameter $a$. Even though the deviation from general relativity for nonzero
values of $a$ can be large, they are still comparable with the variations due
to different modern realistic equations of state. That is why the current
observations of the neutron star masses and radii alone can not put constraints
on the value of the parameter $a$. We also compare our results with those
obtained within the perturbative method and we discuss the differences between
them.
