There are basically, say, 2 x 2 = four forms of results one can extract from NDSolve
, aside from some special-use cases. I think one of these four tends to be more useful than the others, although occasionally I find the other forms useful in particular cases. The one I like is based on a central idiom of Mathematica, the replacement Rule
, and works well in Mathematica. The idea is not common in other languages, and so some may tend to avoid it. The one I like is also in terms of functions $f$ instead of an expression $f(x)$. Despite that everything is an expression in Mathematica, since we're doing calculus here, dealing in terms of functions plays to one of Mathematica's strengths.
The four forms can be viewed as the combinations of two binary divisions. The first division is the choice of the result as a rule or an expression:
NDSolve[eqns, u, {x, xmin, xmax}] (* replacement rule *)
NDSolveValue[eqns, u, {x, xmin, xmax}] (* function expression *)
(* equiv. to u /. NDSolve[eqns, u, {x, xmin, xmax}] *)
The second is the form of u
, either a function or the function-value at the dependent variable x
; for instance
NDSolve[eqns, y, {x, xmin, xmax}] (* (interpolating) function *)
NDSolve[eqns, y[x], {x, xmin, xmax}] (* function-value at x *)
The one I like is
NDSolve[eqns, y, {x, xmin, xmax}]
which returns a rule of the form y -> InterpolatingFunction[<>]
(inside lists, of course). This can be used to inject a function (the solution) into any expression. Mathematica will compute the rest, including any derivatives and integrals of y
that appear in the expression. This is the principal reason I like it.
By contrast
NDSolve[eqns, y[x], {x, xmin, xmax}]
returns a rule of the form y[x] -> InterpolatingFunction[<>][x]
. If it is applied to an expression, it will only replace y[x]
, literally. It will not replace y[2.5]
or y'[x]
, and so forth. This is where the OP ran into trouble with his or her code.
Example
I'll use a slight modification of the OP's example. I'll reduce the number of functions to {z0, z1}
and use slightly more interesting DEs.
Often I find it is convenient to collect all the functions in the differential equation in a list such as vars
below. If vars = {z0, z1}
, then Through[vars[s]]
equals {z0[s], z1[s]}
. The solution solIF
is in the form I prefer, {z0 -> InterpolationFunction[<>],...}
and solS
is in the form {z0 -> InterpolationFunction[<>][s],...}
.
In another side note, the assignment {solIF} =...
works if eqs
are quasilinear system, because the solution will consist of a single set of functions for the dependent variables. If there is an error and NDSolve
does not return a solution, the assignment will fail. One advantage over using solIF = NDSolve[...]
is in code that uses solIF
immediately: On the one hand, if NDSolve
fails, solIF
will equal the original NDSolve[]
call, and every time solIF
is evaluated, it leads to another failed evaluation of NDSolve[]
. On the other hand, if the system is not quasilinear (and I expected it to be so) and NDSolve
returns more than one set of solutions, then the assignment will fail. This lets me knowing I wasn't thinking straight or that the system was not constructed properly. (It's a minor point, but I thought I should explain why I do it.)
(* damped/harmonic oscillators *)
eqs = {z0''[s] + z0'[s]/8 + z0[s] == 0, z1''[s] + 4 z1[s] == 0};
ics = {z0[0] == 0, z0'[0] == 1, z1[0] == 0, z1'[0] == -1};
vars = {z0, z1};
{solIF} = NDSolve[{eqs, ics}, vars, {s, 0, 20}] (* z0 -> ifn, etc *)
{solS} = NDSolve[{eqs, ics}, Through[vars[s]], {s, 0, 20}] (* z0 -> ifn[s], etc *)
To evaluate z0
at 2.5
, say, here's how to use each form. It is somewhat more straightforward with the z0 -> InterpolatingFunction[<>]
form in solIF
.
z0[2.5] /. solIF
z0[s] /. solS /. s -> 2.5
Computing the derivative z1'[s]
is even more natural with solIF
. In the function-value expression form (like the OP's), only the last form of the derivative of z1
is computed correctly, since it is the only code in which the replacement rules are applied to z1[s]
.
z1'[s] /. solS
D[z1[s], s] /. solS
D[z1[s] /. solS, s]
(*
Derivative[1][z1][s]
Derivative[1][z1][s]
InterpolatingFunction[{{0., 20.}}, <>][s]
*)
By contrast, all three form work with the function form:
z1'[s] /. solIF
D[z1[s], s] /. solIF
D[z1[s] /. solIF, s]
Often, we use the solutions to perform other calculations. For instance, we can examine the energy in the oscillator systems represented by {eqs, ics}
. Here is some code for computing the energy:
ClearAll[energy];
SetAttributes[energy, Listable];
energy[harmonicODE_, s0_] := Module[{z, s, m, γ, k},
{{z, s}} = Cases[harmonicODE, Derivative[2][z_][s_] :> {z, s}, Infinity, 1];
{k, γ, m} = Coefficient[harmonicODE /. Equal -> Subtract, {z[s], z'[s], z''[s]}];
k z[s0]^2/2 + m z'[s0]^2/2
];
energy[eqs, s]
(*
{z0[s]^2/2 + 1/2 z0'[s]^2,
2 z1[s]^2 + 1/2 z1'[s]^2}
*)
These are symbolic expressions for the energies. We can evaluate them by injecting the solution functions in solIF
into the energies and plugging in values for s
, say, via Plot[]
:
Plot[energy[eqs, s] /. solIF // Evaluate, {s, 0, 20}]

I don't know a simple way to do that with solS
.
I won't say much about NDSolveValue[]
. It returns InterpolatingFunction[<>]
or InterpolatingFunction[<>][s]
. It can be convenient if you just want to evaluate the function. One can inject them in an expression by writing your own rules, but that can be a bit inconvenient if you have several variable {z0, z1,...}
as in the OP's case. You're probably better off using NDSolve[]
in that case.