  We revisit a method to incorporate the Vainshtein screening mechanism in
N-body simulations proposed by R. Scoccimarro in~\cite{Scoccimarro:2009eu}. We
further extend this method to cover a subset of Horndeski theories that evade
the bound on the speed of gravitational waves set by the binary neutron star
merger GW170817. The procedure consists of the computation of an effective
gravitational coupling that is time and scale dependent, $G_{\rm
eff}\left(k,z\right)$, where the scale dependence will incorporate the
screening of the fifth-force. This is a fast procedure that when contrasted to
the alternative of solving the full equation of motion for the scalar field
inside N-body codes, reduces considerably the computational time and complexity
required to run simulations. To test the validity of this approach in the
non-linear regime, we have implemented it in a COmoving Lagrangian
Approximation (COLA) N-body code, and ran simulations for two gravity models
that have full N-body simulation outputs available in the literature, nDGP and
Cubic Galileon. We validate the combination of the COLA method with this
implementation of the Vainshtein mechanism with full N-body simulations for
predicting the boost function: the ratio between the modified gravity
non-linear matter power spectrum and its General Relativity counterpart. This
quantity is of great importance for building emulators in beyond-$\Lambda$CDM
models, and we find that the method described in this work has an agreement of
below $2\%$ for scales down to $k \approx 3h/$Mpc with respect to full N-body
simulations.
