In this talk I discuss recent joint work with Oksana Chkrebtii, Prof. Dave Campbell and Prof. Mark Girolami, in which we develop a probabilistic formalism for solving systems of differential equations. This generalises classes of existing numerical solvers while making the modelling assumptions explicit. The approach discussed yields a probability distribution on a function space of possible solutions, instead of a single deterministic solution that approximately satisfies model dynamics to within a given error tolerance. Viewing solution estimation as an inference problem (O'Hagan, 1992; Skilling, 1991) allows us to quantify solver error using the tools of Bayesian function estimation. In particular, we make use of Gaussian process priors on an underlying function space, while incorporating regularity assumptions by modelling states and their derivatives by their kernel integral transforms.