  In this work by using a numerical analysis, we investigate in a quantitative
way the late-time dynamics of scalar coupled $f(R,\mathcal{G})$ gravity.
Particularly, we consider a Gauss-Bonnet term coupled to the scalar field
coupling function $\xi(\phi)$, and we study three types of models, one with
$f(R)$ terms that are known to provide a viable late-time phenomenology, and
two Einstein-Gauss-Bonnet types of models. Our aim is to write the Friedmann
equation in terms of appropriate statefinder quantities frequently used in the
literature, and we numerically solve it by using physically motivated initial
conditions. In the case that $f(R)$ gravity terms are present, the contribution
of the Gauss-Bonnet related terms is minor, as we actually expected. This
result is robust against changes in the initial conditions of the scalar field,
and the reason is the dominating parts of the $f(R)$ gravity sector at late
times. In the Einstein-Gauss-Bonnet type of models, we examine two distinct
scenarios, firstly by choosing freely the scalar potential and the scalar
Gauss-Bonnet coupling $\xi(\phi)$, in which case the resulting phenomenology is
compatible with the latest Planck data and mimics the
$\Lambda$-Cold-Dark-Matter model. In the second case, since there is no
fundamental particle physics reason for the graviton to change its mass, we
assume that primordially the tensor perturbations propagate with the speed
equal to that of light's, and thus this constraint restricts the functional
form of the scalar coupling function $\xi(\phi)$, which must satisfy the
differential equation $\ddot{\xi}=H\dot{\xi}$.
