We develop a full four-dimensional numerical code to study scalar gravitational radiation emitted from binary systems and probe the Vainshtein mechanism in situations that break the static and spherical symmetry, relevant for binary pulsars as well as black holes and neutron stars binaries. The present study focuses on the cubic Galileon which arises as the decoupling limit of massive theories of gravity. Limitations associated with the numerical methods prevent us from reaching a physically realistic hierarchy of scales; nevertheless, within this context we observe the same power law scaling of the radiated power as previous analytic estimates, and confirm a strong suppression of the power emitted in the monopole and dipole as compared with quadrupole radiation. Following the trend to more physically realistic parameters, we confirm the suppression of the power emitted in scalar gravitational radiation and the recovery of general relativity with good accuracy. This paves the way for future numerical work, probing more generic, physically relevant situations and sets of interactions that may exhibit the Vainshtein mechanism.