# Numerically solving ODEs — how to estimate the solution between the nodes?

I have heard about a lot of fancy numerical methods for solving ODEs. I know that there are methods that (assuming sufficient smoothness) asymptotically give a low error, like the Runge-Kutta methods. These estimate the solution in a set of points $$t_0$$, $$t_1$$, etc. But what if I want to have a function that is close to the correct solution everywhere, not just in a discrete set of points?

I can extend the numerical solution to a piecewise linear function. This will be a continuous function and it will converge to the correct solution if the step-size goes to zero.

But the error estimate will be poor in most places unless I use a very small step-size, which rather defeats the purpose of using a high-order method. So how does one go about estimating the solution in practice between the $$t_i$$?

Remember that whatever method you use for solving $y^{\prime}=f(t,y)$, be it Runge-Kutta, Bulirsch-Stoer (extrapolative), Gear/Adams multistep, or fancier methods, one always has a triple of values $(t_i,y_i,y_i^{\prime})$ ($y_i^{\prime}=f(t_i,y_i)$) at each evaluation point. Thus, one can always do cubic Hermite interpolation across the points $(t_i,y_i)$ and $(t_{i+1},y_{i+1})$ (note that I am not assuming that the evaluation points are equispaced, as is often the case when doing adaptive solving). If the underlying method is at most third order accurate, cubic Hermite is a good choice.
Now, the upshot is that modern implementations of DE solvers always support so-called "dense output"; briefly, associated with a method with $p$th order "accuracy" (in quotes since "high order doesn't always imply high accuracy" ;) ) has with it an associated $p$th order interpolating function. To use my favorite example, the $(4,5)$ adaptive Runge-Kutta solver based on the Dormand-Prince coefficients has with it an associated fifth order interpolating function. The special properties inherent in the coefficients allow the existence of an associated interpolating function with the same order of "accuracy"; in general, not all Runge-Kutta coefficients will have an associated "nice" interpolating function (but again, one can always do cubic Hermite).
Take for example Runge-Kutta methods. Let $t_1$, $t_2$, etc., be the points given by the step-size of the method. Then the Runge-Kutta method not only give you the approximate solution at points $t_1$, $t_2$, etc., but also at intermediate points between, say, $t_1$ and $t_2$, depending on the order of the method. You can use these intermediate points (between $t_1$ and $t_2$, together with the endpoints $t_1$ and $t_2$) to construct an interpolation polynomial between $t_1$ and $t_2$, so that the resulting approximate solution will be a piecewise polynomial function, more general than piecewise linears, and the accuracy will as good as the accuracy of the Runge-Kutta method at the mesh points.