We study the late cosmological evolution, from the nonrelativistic matter dominated era to the dark energy era, in modified gravity models described by Degenerate Higher-Order Scalar-Tensor (DHOST) theories. They represent the most general scalar-tensor theories propagating a single scalar degree of freedom and include Horndeski and Beyond Horndeski theories. We provide the homogeneous evolution equations for any quadratic DHOST theory, without restricting ourselves to theories where the speed of gravitational waves coincides with that of light since the present constraints apply to wavelengths much smaller than cosmological scales. To illustrate the potential richness of the cosmological background evolution in these theories, we consider a simple family of shift-symmetric models, characterized by three parameters and compute the evolution of dark energy and of its equation of state. We also identify the regions in parameter space where the models are perturbatively stable.

We study a dark energy cubic bi-Galileons field model based on truncation of the recently proposed generalized covariant multi-Galileons model. We investigate the cosmological dynamic of the model by the theory of dynamical systems through the analysis of the properties of the fixed points in each cosmological epoch. We show the existence of two tracker solutions, one of which is that of the cubic single Galileon model and the other solution is the signature of the second Galileon field. Exploiting the competition between the two Galileon fields, we find a dark energy solution that avoids the approach to the tracker solution with dark energy equation of state $w_{DE}=-2$ during the matter epoch which is disfavored by the observational data. We study also the growth rate of matter perturbations. Using recent $f\sigma_{8}$ redshift space distortion (RSD) and model-independent observational Hubble (OHD) data sets, we place observational constraints on the coupling constant and cosmological parameters of the bi-Galileons model through Monte Carlo numerical method based on the Metropolis-Hastings algorithm. We find that the amplitude of growth matter fluctuations is consistent with the Planck15 data and ease the tension between early and later clustering, and fits better the data from the DES survey over the data from KiDS-450 survey. We also find that the best fit value for the Hubble constant is compatible with new measurements of Cepheid-supernovae distance scale. Finally, we perform a model selection through the Bayes factor and found that the bi-Galileons model is disfavored in comparison to the $\Lambda$CDM model, but slightly preferred to $w$CDM model.