  This paper analyses the cosmological consequences of a modified theory of
gravity whose action integral is built from a linear combination of the Ricci
scalar $R$ and a quadratic term in the covariant derivative of $R$. The
resulting Friedmann equations are of the fifth-order in the Hubble function.
These equations are solved numerically for a flat space section geometry and
pressureless matter. The cosmological parameters of the higher-order model are
fit using SN Ia data and X-ray gas mass fraction in galaxy clusters. The
best-fit present-day $t_{0}$ values for the deceleration parameter, jerk and
snap are given. The coupling constant $\beta$\ of the model is not univocally
determined by the data fit, but partially constrained by it. Density parameter
$\Omega_{m0}$ is also determined and shows weak correlation with the other
parameters. The model allows for two possible future scenarios: there may be
either a premature Big Rip or a Rebouncing event depending on the set of values
in the space of parameters. The analysis towards the past performed with the
best-fit parameters shows that the model is not able to accommodate a
matter-dominated stage required to the formation of structure.
