  We present and test a new numerical method to determine second-order
Lagrangian displacement fields in the context of modified gravity (MG)
theories. We start from the extension of Lagrangian Perturbation Theory to a
class of MG models that can be described by a parametrized Poisson equation,
with the introduction of a scale-dependent function. We exploit fast Fourier
transforms to compute the full source term of the differential equation for the
second-order Lagrangian displacement field. We compare its mean to the source
term computed for specific configurations for which a k-dependent solution can
be found numerically. We choose the configuration that best matches the full
source term, thus obtaining an approximate factorization of the second-order
displacement field as the space term valid for standard gravity times a
k-dependent, second-order growth factor $D_2(k,t)$. This approximation is used
to compute second order displacements for particles. The method is tested
against N-body simulations run with standard and $f(R)$ gravity: we rely on the
results of a friends-of-friends code run on the N-body snapshots to assign
particles to halos, then compute the halo power spectrum. We find very
consistent results for the two gravity theories: second-order LPT (2LPT) allows
to recover the halo power spectrum of N-body simulations within $\sim 10\%$
precision to $k\sim 0.2-0.4\ h\ {\rm Mpc}^{-1}$, as well as halo positions,
with an error that is a fraction of the inter-particle distance. We show that,
when considering the same level of non-linearity in the density field, the
performance of 2LPT with MG is the same (within $1\%$) as the one obtained for
the standard $\Lambda$CDM model with General Relativity. When implemented in a
computer code, this formulation of 2LPT can quickly generate dark matter
distributions with $f(R)$ gravity, and can easily be extended to other MG
theories.
