  We study the structure of neutron stars in perturbative f(R) gravity models
with realistic equations of state. We obtain mass-radius relations in a gravity
model of the form f(R)=R+\alpha R^2. We find that deviations from the results
of general relativity, comparable to the variations due to using different
equations of state (EoS'), are induced for |alpha| ~ 10^9 cm^2. Some of the
soft EoS' that are excluded within the framework of general relativity can be
reconciled with the 2 solar mass neutron star recently observed for certain
values of alpha within this range. For some of the EoS' we find that a new
solution branch, which allows highly massive neutron stars, exists for values
of alpha greater than a few 10^9 cm^2. We find constraints on alpha for a
variety of EoS' using the recent observational constraints on the mass-radius
relation. These are all 5 orders of magnitude smaller than the recent
constraint obtained via Gravity Probe B for this gravity model. The associated
length scale \sqrt{alpha} ~ 10^5 cm is only an order of magnitude smaller than
the typical radius of a neutron star, the probe used in this test. This implies
that real deviations from general relativity can be even smaller.
