We describe a set of Gaussian Process based approaches that can be used to
solve non-linear Ordinary Differential Equations. We suggest an explicit
probabilistic solver and two implicit methods, one analogous to Picard
iteration and the other to gradient matching. All methods have greater accuracy
than previously suggested Gaussian Process approaches. We also suggest a
general approach that can yield error estimates from any standard ODE solver