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+αR². 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 |α| ∼ 10⁹ cm². 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 α 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 α greater than a few 10⁹ cm². We find constraints on α 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 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.