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

Question

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. >>

How can I solve the problem? Please help.

Answer 1

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.