Résumé : 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.