ABSTRACT We present MG-evolution, an N-body code simulating the cosmological structure formation for parametrized modifications of gravity. It is built from the combination of parametrized linear theory with a parametrization of the deeply non-linear cosmological regime extrapolated from modified spherical collapse computations that cover the range of known screening mechanisms. We test MG-evolution, which runs at the speed of conventional ΛCDM simulations, against a suit of existing exact model-specific codes, encompassing linearized and chameleon f(R) gravity as well as the normal branch of the Dvali–Gabadadz–Porrati braneworld model, hence covering both large-field value and large-derivative screening effects. We compare the non-linear power spectra produced by the parametrized and model-specific approaches over the full range of scales set by the box size and resolution of our simulations, k = (0.05 − 2.5) $h\, \mathrm{Mpc}^{-1}$, and for two redshift slices, z = 0 and z = 1. We find sub-percent to one-percent level recovery of all the power spectra generated with the model-specific codes for the full range of scales. MG-evolution can be used for generalized and accurate tests of gravity and dark energy with the increasing wealth of high-precision cosmological survey data becoming available over the next decade.