First, InterpolationOrder can be set to integer and make a difference, but there are limitations. Whether the settings other than Automatic and All are useful is another question. For the default Method -> Automatic on an ordinary ODE, which uses LSODA, these are the only options and other settings revert to one of these. Perhaps the main thing to remember is that InterpolationOrder works in conjunction with the integration method of NDSolve and affects the kind of solution constructed; it does not have the same function as it does in Interpolation.
Numerical integration of differential equations might be thought of as consisting of two parts, the first part being the construction of discrete steps of the solution and the other part using interpolation to fill in the gaps. The accuracy of the first part is controlled by PrecisionGoal and AccuracyGoal. Let's consider why one might want interpolation orders above, below or in between the orders implied by Automatic and All.
The basic purpose of raising the interpolation order is to reduce the interpolation error between the steps taken by NDSolve. A higher order interpolation requires more memory to store the result and more time to evaluate the interpolation. It is normally impossible to improve the accuracy beyond the underlying order of the integration method of NDSolve, but the see examples below for an exception. Since most integration orders are low, it is unclear that there will often be an advantage in requesting an interpolation order between the Automatic and All. I can imagine cases in which the user would be satisfied to save memory and get just the steps, i.e., a linear interpolation; however, NDSolve does not support this. I do not see any other advantage in lowering the order below 3. Note that since the default method produces an accurate-looking graph, raising the order usually will not result in a visible difference. It should be used when one wants to calculate accurately with the solution, especially when using its derivatives as the examples below show.
It seems that if["InterpolationOrder"] does not always give reliable information for an InterpolatingFunction if. We cannot rely on it to test whether NDSolve is returning a solution of the requested interpolation order. The integration methods used by NDSolve may have restrictions on both how small or how great the interpolation order can be. There are broadly two kinds of InterpolatingFunction that NDSolve might return. One uses piecewise Hermite interpolation, in which the interpolation order is determined automatically by the local data (function and derivative values) at each step. The other uses piecewise series. Of the series, it generally uses either local Taylor series or Chebyshev series; for small-order steps, it will use cubic Hermite interpolation.
For instance, for an ODE of order $n$, it seems wasteful not to use at a minimum Hermite interpolation of order $2n+1$. This is because at each step, NDSolve calculates the $n+1$ values of
$$y_k, y'_k, y''_k, \dots, y^{(n)}_k$$
and it might as well store them in the answer, unless the memory requirements are prohibitive. This also means the solution will have the same order of smoothness as the ODE. So $2n+1$ seems to be a lower bound on the interpolation order, since the $2\times(n+1)$ values at both end points of each step determine polynomial of degree $2n+1$. (In fact, NDSolve seems to do this only up to the first 12 derivatives.) For first-order ODEs, the lower bound is 3. If it is desired to have a lower-order, worse solution, then re-interpolate; the InterpolatingFunction tools make it easy, but I'm not sure why you would do that. (However, beware the easy way will unpack the packed interpolation data. It is difficult to do this without unpacking it.) Many integration methods switch to local Taylor series or Chebyshev series for piecewise interpolation when the requested interpolation order is higher (than 3 most commonly or perhaps always). The ones that use Chebyshev series seem to honor arbitrarily high requests, even though it usually does not improve quality once the interpolation order is higher than the underlying order of the method.
Remarks on some methods:
"LSODA" (the default or Automatic method): This method has just two interpolation options, which perhaps makes the option seem nearly useless. The default Hermite (when InterpolationOrder is 3 or less or Automatic) and local Taylor series. It uses series of varying orders, sometimes less than the default Hermite method for ODEs of order two and higher.
"ExplicitRungeKutta", "ImplicitRungeKutta", "Extrapolation", "DoubleStep": These use Chebyshev series. InterpolationOrder -> All yields varying orders for "Extrapolation"; for the others, with the default "DifferenceOrder", the interpolation is equivalent to order 9.
"ExplicitEuler": This has just two options, the default Hermite, or, when the order is set higher than 3, cubic Hermite even if the default Hermite has a higher order. Go figure. Used as a submethod of "DoubleStep", "ExplicitEuler" yields Chebyshev series of the requested order. Judging from the coefficients, the true order is 2, no matter what the requested order.
Other methods will show a similar range of behaviors.
Some examples:
First some utilities. We call NDSolve often and repeatedly on the same. Memoizing is a simple way to compute each unique call only once.
mem : ndsolve[stuff___] := (* memoize NDSolve calls *)
With[{res = NDSolve[stuff]},
If[FreeQ[res, NDSolve],
mem = First@res, (* assumes there is just one solution *)
First@res]
];
(* OP's first-order 2D example *)
ode1 = {x'[t] == y[t], y'[t] == -x[t], x[0] == 1, y[0] == 0};
sol1[Infinity] = First@DSolve[ode1, {x, y}, t];(* exact solution {Cos[t], -Sin[t]} *)
sol1[n_, meth_] :=
ndsolve[ode1, {x, y}, {t, 0, 10}, InterpolationOrder -> n, Method -> meth];
(* equivalent second-order example *)
ode2 = {x''[t] == -x[t], x[0] == 1, x'[0] == 0};
sol2[Infinity] = First@DSolve[ode2, x, t];(* exact solution Cos[t] *)
sol2[n_, meth_] :=
ndsolve[ode2, x, {t, 0, 10}, InterpolationOrder -> n, Method -> meth];
irkmeth = {"ImplicitRungeKutta", (* "ImplicitRungeKutta" method used below *)
"Coefficients" -> NDSolve`ImplicitRungeKuttaGaussCoefficients,
DifferenceOrder -> 9};
Results for various methods on ode2 with different InterpolationOrder settings. The O[n] indicates order n. The coefficient indicates how many steps have that order of the corresponding type.
The following shows how the solution converges toward the exact solution as InterpolationOrder increases in a Chebyshev-series interpolation.
The exponent (base 10) of the error of an implicit Runge-Kutta method with Gauss coefficients of order 9 and machine precision; the interpolation order is shown in the SetterBar in the label. One can see that InterpolationOrder -> 4 has worse error than InterpolationOrder -> Automatic except for its higher derivatives derivatives. One can see that InterpolationOrder -> 5 has about the same error for x[t] but better error for its derivatives.
With higher precision, something interesting happens in this example. When the interpolation order exceeds All (or 9), the quality of the derivatives continues to improve up to order 11.
The exponent (base 10) of the error of an implicit Runge-Kutta method with Gauss coefficients of order 9 and working precision 32; the interpolation order is shown in the SetterBar in the label. One can see that the accuracy of the derivatives improves when the interpolation order increases up to order 11, past the order 9 that is equivalent to All.
