  We describe and demonstrate the potential of a new and very efficient method
for simulating certain classes of modified gravity theories, such as the widely
studied $f(R)$ gravity models. High resolution simulations for such models are
currently very slow due to the highly nonlinear partial differential equation
that needs to be solved exactly to predict the modified gravitational force.
This nonlinearity is partly inherent, but is also exacerbated by the specific
numerical algorithm used, which employs a variable redefinition to prevent
numerical instabilities. The standard Newton-Gauss-Seidel iterative method used
to tackle this problem has a poor convergence rate. Our new method not only
avoids this, but also allows the discretised equation to be written in a form
that is analytically solvable. We show that this new method greatly improves
the performance and efficiency of $f(R)$ simulations. For example, a test
simulation with $512^3$ particles in a box of size $512 \, \mathrm{Mpc}/h$ is
now 5 times faster than before, while a Millennium-resolution simulation for
$f(R)$ gravity is estimated to be more than 20 times faster than with the old
method. Our new implementation will be particularly useful for running very
high resolution, large-sized simulations which, to date, are only possible for
the standard model, and also makes it feasible to run large numbers of lower
resolution simulations for covariance analyses. We hope that the method will
bring us to a new era for precision cosmological tests of gravity.
