  The interpretation of numerical methods, such as finite difference methods
for differential equations, as point estimators allows for formal statistical
quantification of the error due to discretisation in the numerical context.
Competing statistical paradigms can be considered and Bayesian probabilistic
numerical methods (PNMs) are obtained when Bayesian statistical principles are
deployed. Bayesian PNM are closed under composition, such that uncertainty due
to different sources of discretisation can be jointly modelled and rigorously
propagated. However, we argue that no strictly Bayesian PNM for the numerical
solution of ordinary differential equations (ODEs) have yet been developed. To
address this gap, we work at a foundational level, where a novel Bayesian PNM
is proposed as a proof-of-concept. Our proposal is a synthesis of classical Lie
group methods, to exploit the underlying structure of the gradient field, and
non-parametric regression in a transformed solution space for the ODE. The
procedure is presented in detail for first order ODEs and relies on a certain
technical condition -- existence of a solvable Lie algebra -- being satisfied.
Numerical illustrations are provided.
