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.