In this talk I will outline a fast and accurate approach to introduce modified gravity (MG) consistently in Cosmological N-Body simulations, using the COLA (COmoving Lagrangian Acceleration) method, a fast and approximate method to generate non-linear realizations of the density field. The matter power spectrum computed from these simulations is shown to be in agreement with known results in the literature at all scales, i.e., from the largest scales, where relativistic effects dominate the signal, to the smallest scales, where screening is already active and shielding the fifth force introduced by the scalar field.