  We introduce and demonstrate the power of a method to speed up current
iterative techniques for N-body modified gravity simulations. Our method is
based on the observation that the accuracy of the final result is not
compromised if the calculation of the fifth force becomes less accurate, but
substantially faster, in high-density regions where it is weak due to
screening. We focus on the nDGP model which employs Vainshtein screening, and
test our method by running AMR simulations in which the solutions on the finer
levels of the mesh (high density) are not obtained iteratively, but instead
interpolated from coarser levels. We show that the impact this has on the
matter power spectrum is below $1\%$ for $k < 5h/{\rm Mpc}$ at $z = 0$, and
even smaller at higher redshift. The impact on halo properties is also small
($\lesssim 3\%$ for abundance, profiles, mass; and $\lesssim 0.05\%$ for
positions and velocities). The method can boost the performance of modified
gravity simulations by more than a factor of 10, which allows them to be pushed
to resolution levels that were previously hard to achieve.
