  We present numerical solutions of several spacetimes of physical interest,
including binary black hole mergers, in shift-symmetric
Einstein-scalar-Gauss-Bonnet (ESGB) gravity, and describe our methods for
solving the full equations of motion, without approximation, for general
spacetimes. While we concentrate on the specific example of shift-symmetric
ESGB, our methods, which make use of a recently proposed modification to the
generalized harmonic formulation, should be generally applicable to all
Horndeski theories of gravity (including general relativity). We demonstrate
that these methods can stably follow the formation of scalar clouds about
initially vacuum non-spinning and spinning black holes for values of the
Gauss-Bonnet coupling approaching the maximum value above which the
hyperbolicity of the theory breaks down in spherical symmetry. We study the
collision of black holes with scalar hair, finding that the theory remains
hyperbolic in the spacetime region exterior to the black hole horizons in a
similar regime, which includes cases where the deviations from general
relativity in the gravitational radiation is appreciable. Finally, we
demonstrate that these methods can be used to follow the inspiral and merger of
binary black holes in full ESGB gravity. This allows for making predictions for
Horndeski theories of gravity in the strong-field and non-perturbative regime,
which can confronted with gravitational wave observations, and compared to
approximate treatments of modifications to general relativity.
