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

You can then use your favorite interpolation technique. As you say, you have a bunch of points (t,x). So you can feed them to a spline routine, or a polynomial fit, or whatever. If you have some knowledge of the functional form, you can take these data points to fit the form, but it may be better to use your differential equation for that.

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.

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).

I could say more, but the books by Hairer/Nørsett/Wanner have an extensive discussion on dense output (and they say it better than what I can hope to say), as well as usable routines (also available within the site). You would do well to study them.