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.