Using NDSolve initial conditions as parameters for a function and then a replacement rule

Question

I am trying to make a function based on two initial conditions of a set of differential equations:

sol[a, b] := 
NDSolve[{x1'[t] == ..., y1'[t] == ..., z1'[t] == ..., 
         x2'[t] == ..., y2'[t] == ..., z2'[t] == ...,
         p1'[t] == ..., p2'[t] == ..., p3'[t] == ...,
         e1'[t] == ..., e2'[t] == ..., e3'[t] == ...
         x2[0] == a, y2[0] == b (*insert other initial conditions here*)},
{x2[t], y2[t], z2[t], x1[t], y1[t], z1[t], p1[t], p2[t], p3[t], 
e1[t], e2[t], e3[t]}, {t, 0, T}];

However, this gives the error "NDSolve::ndinnt: Initial condition a is not a number or a rectangular array of numbers."

I also want to use this function as a replacement rule in order to plot a contour graph giving the value of e1[c] (for a certain fixed c) based on variables a and b:

(*c is some constant*)
g = e1[c] /. sol2[a, b];
ContourPlot[g, {a, -1, 1}, {b, -1, 1}]

But, this results in another error "ReplaceAll::reps: NDSolve[...] is neither a list of replacement rules nor a valid dispatch table, and so cannot be used for replacing."

Any way to explain (and hopefully fix) these two errors?

Answer 1

As mentioned by bbgodfrey in a comment, you can use ParametricNDSolve. As a toy example, let's do the following:

sols = ParametricNDSolve[{
  x'[t] == -x[t]
  , y'[t] == -y[t]
  , e'[t] == -e[t] + x[t]/y[t]
  , x[0] == a, y[0] == b , e[0] == 0
 }
 , {x, y, e}
 , {t, 0, 5}
 , {a, b}
]

Then

ContourPlot[Evaluate[e[a, b][c] /. sols /. c -> 1], {a, 0, 1}, {b, 0, 1}]

results in

Notes about OP's code

You can actually make your method work. It's not advisable, since it will be much slower due to calling NDSolve for every unique pair of a and b. Nonetheless, here's how you would go about doing it (included for the purpose of learning more about MMA by explaining what went wrong):

Clear[sols]
sols[a_, b_] := NDSolve[{x'[t] == -x[t]
  , y'[t] == -y[t]
  , e'[t] == -e[t] + x[t]/y[t]
  , x[0] == a, y[0] == b , e[0] == 0
 }
 , {x, y, e}
 , {t, 0, 5}
];
ContourPlot[e[1] /. sols[a, b], {a, 0, 1}, {b, 0, 1}]

We define the function using SetDelayed (:=) due to the fact that NDSolve requires values for a and b in order to run. By using Set (=), MMA will try to run NDSolve immediately, even though a and b aren't defined as numbers. This is the reason for the first error.

We have also included underscores in the definition of sols: sols[a_, b_] instead of sols[a, b]. This is the way that MMA "functions" are defined. If you use a and b in the definition, then the right-hand side will be set to sols[a, b] literally, and you will not be able to change the values of a and b. Using the Patterns a_ and b_ tells MMA that whatever you put in place of a_ and b_ will be fed to a and b in the right-hand-side expression.

As an example, look at the output when you run the following:

Clear[f, a, b]
f[a, b] = a + b
f[1, 2]
Clear[f, a, b]
f[a_, b_] = a + b
f[1, 2]

Finally, the second error comes about because you tried to evaluate NDSolve without specifying values for a and b. When you do this, MMA will just return the original expression and set it equal to sols[a,b]. Thus, when you try to use it as a replacement rule in your call to ContourPlot, it's not understood because expressions with NDSolve as a Head are not replacement rules.