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.