# 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] == 0, z[0] == 0, y'[0] == -0.1, z'[0] == 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]==0,z[0]==0,(y^[Prime])[0]==-0.1,(z^[Prime])[0]==0.05},{y,z},{x,0,10}]] are not the same shape. >>

• What's wrong with ParametricNDSolveValue? Oct 16, 2017 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". Oct 16, 2017 at 13:10
• I am not sure about the // issue, using First[] also yields the same error. Oct 16, 2017 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. Oct 16, 2017 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. Oct 16, 2017 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] == 0, z[0] == 0, y'[0] == -0.1, z'[0] == 0.05}, {y, z}, {x,
0, 10}]

funs = sol[1]

Through[funs[1]]
(* {-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] == 0, z[0] == 0, y'[0] == -0.1,
z'[0] == 0.05}, {y, z}, {x, 0, 10}, a]


Now psol[1] returns the same thing as sol[1].

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. Oct 16, 2017 at 15:39