Why doesn't ?NumericQ work in NDSolve with two variable?

I am trying to solve the following set of differential eqiations

$\qquad y''(x) + a\,\sin(x) z(x) = 0 \quad z''(x) + a\,z(x) y(x) = 0$

where $a$ is a parameter. I want to get a solution for many values ofd $a$, so I constructed something like the following (I didn't want to use ParametricNDSolve as the actual equation I need to solve is stiff as well as I need to solve some algebraic equation in terms of solutions to find out a specific $a$)

{y1[a_?NumericQ], z1[a_?NumericQ]} :=
{y1[a],
z1[a]} = {y, z} /.
NDSolve[
{y''[x] + a*Sin[x] z[x] == 0, z''[x] + a*z[x]
y[x] == 0, y == 0, z == 0, y' == -0.1, z' == 0.05},
{y, z}, {x, 0, 10}] // First

I have used these type of functionalized constructs for single differential equation, but in this case of system of equations the following error is showing

SetDelayed::shape: Lists {y1[a_?NumericQ],z1[a_?NumericQ]} and {y1[a],z1[a]}=First[{y,z}/. NDSolve[{(<<1>>^(<<1>>))[<<1>>]+Times[<<3>>]==0,(<<1>>^(<<1>>))[<<1>>]+Times[<<3>>]==0,y==0,z==0,(y^[Prime])==-0.1,(z^[Prime])==0.05},{y,z},{x,0,10}]] are not the same shape. >>

• What's wrong with ParametricNDSolveValue? – Szabolcs Oct 16 '17 at 13:03
• The original equation which I want to solve is of three variables besides it is a stiff equation, i.e., I need to use Method -> "StiffnessSwitching". – Archimedes Oct 16 '17 at 13:10
• I am not sure about the // issue, using First[] also yields the same error. – Archimedes Oct 16 '17 at 13:12
• I was wrong about First. The problem is that there are in fact two separate definitions, and they must be independent. I would have a single function that returns two values, or similar, if I really wanted memoization. – Szabolcs Oct 16 '17 at 13:18
• Does ParametricNDSolveValue not work with "StiffnessSwitching"? I just tried it and it worked fine. Furthermore, it seems to support solution caching too, which means that you do not need memoization. – Szabolcs Oct 16 '17 at 13:19

You have an expression of the form

{a, b} := Set[...]

Assigning multiple values simultaneously requires lists of the same shape to be explicitly present on both sides of the :=, i.e. {a,b} = {c,d}. This is not the case here, hence the error.

You can't have an explicit list on the RHS because you use a single computation to get both functions, not two computations.

If you want to avoid ParametricNDSolveValue, define a single function instead, which returns two interpolating functions.

sol[a_?NumericQ] :=
NDSolveValue[{y''[x] + a*Sin[x] z[x] == 0, z''[x] + a*z[x] y[x] == 0,
y == 0, z == 0, y' == -0.1, z' == 0.05}, {y, z}, {x,
0, 10}]

funs = sol

Through[funs]
(* {-0.103905, 0.0504225} *)

That said, I strongly suggest just using ParametricNDSolveValue, which was made exactly for this purpose.

psol = ParametricNDSolveValue[{y''[x] + a*Sin[x] z[x] == 0,
z''[x] + a*z[x] y[x] == 0, y == 0, z == 0, y' == -0.1,
z' == 0.05}, {y, z}, {x, 0, 10}, a]

Now psol returns the same thing as sol.

ParametricNDSolveValue also gives you caching for free, so you do not have to try to implement your own memoization. See the "ParametricCaching" option.

• Thanks for the detailed answer. – Archimedes Oct 16 '17 at 15:39