  In this paper, we treat the problem of evaluating the asymptotic error in a
numerical integration scheme as one with inherent uncertainty. Adding to the
growing field of probabilistic numerics, we show that Gaussian process
regression (GPR) can be embedded into a numerical integration scheme to allow
for (i) robust selection of the adaptive step-size parameter and; (ii)
uncertainty quantification in predictions of putatively converged numerical
solutions. We present two examples of our approach using Richardson's
extrapolation technique and the Bulirsch-Stoer algorithm. In scenarios where
the error-surface is smooth and bounded, our proposed approach can match the
results of the traditional polynomial (parametric) extrapolation methods. In
scenarios where the error surface is not well approximated by a finite-order
polynomial, e.g. in the vicinity of a pole or in the assessment of a chaotic
system, traditional methods can fail, however, the non-parametric GPR approach
demonstrates the potential to continue to furnish reasonable solutions in these
situations.
