We investigate the effect of a strong magnetic field on the structure of neutron stars in a model with perturbative $f(R)$ gravity. The effect of an interior strong magnetic field of about $10^{17 \sim 18}$ G on the equation of state is derived in the context of a quantum hadrodynamics (QHD) model. We solve the modified spherically symmetric hydrostatic equilibrium equations derived for a gravity model with $f(R)=R+\alpha R^2$. Effects of both the finite magnetic field and the modified gravity are detailed for various values of the magnetic field and the perturbation parameter $\alpha$ along with a discussion of their physical implications. We show that there exists a parameter space of the modified gravity and the magnetic field strength, in which even a soft equation of state can accommodate a large ($> 2$ M$_\odot$) maximum neutron star mass through the modified mass-radius relation.

We study the structure of neutron stars in R+\beta\ R^{\mu \nu} R_{\mu \nu} gravity model with perturbative method. We obtain mass--radius relations for six representative equations of state (EoSs). We find that, for |\beta| ~ 10^11 cm^2, the results differ substantially from the results of general relativity. Some of the soft EoSs that are excluded within the framework of general relativity can be reconciled for certain values of \beta\ of this order with the 2 solar mass neutron star recently observed. For values of \beta\ greater than a few 10^11 cm^2 we find a new solution branch allowing highly massive neutron stars. By referring some recent observational constraints on the mass--radius relation we try to constrain the value of \beta\ for each EoS. The associated length scale \sqrt{\beta} ~ 10^6 cm is of the order of the typical radius of neutron stars implying that this is the smallest value we could find by using neutron stars as a probe. We thus conclude that the true value of \beta\ is most likely much smaller than 10^11 cm^2.