Abstract: Differential equation models offer a succinct way of describing rates of change using few but readily interpretable parameters. In most interesting cases analytic solutions do not exist and numerical methods are used to approximate a solution over a discretization grid. We explore the use of probability models for uncertainty arising from the discretization of ordinary or partial differential equation solutions. Viewing the system solution as an inference problem allows us to quantify numerical uncertainty using the tools of smoothing and gaussian process regression. A formalism for inferring differential equation model trajectories is proposed through a Bayesian updating scheme based on interrogations of the model derivatives. The approach provides estimates for the differential equation model solution and derivatives thereof which can be incorporated into the inference problem. The proposed approach is demonstrated to capture the functional structure and magnitude of the discretization error, while attaining computational scaling of the same order as standard numerical solver methods. This allows us to define a trade-off between accuracy and discretization grid size. Our approach is demonstrated on ordinary and partial differential equation models, ill-conditioned mixed boundary value problems, and delay differential equations. This talk represents joint work with Oksana Chkrebtii, Mark Girolami, and Ben Calderhead and appeared as a discussion paper in Bayesian Analysis in December 2016.