# Error when using NDSolveValue. Lists are not the same shape

Please, help me with the following issue. I am solving the following system of the coupled PDEs

with the following boundary conditions

and initial conditions.

After performing the following simplified code (for simplicity I set all parameters equal to 1):

(*specify the system*)
op =
{-Inactive[Grad][n[t, x], t] + Laplacian[n[t, x], x] -
n[t, x] - (n[t, x])^2,
-Inactive[Grad][a[t, x], t] + Inactive[Laplacian][a[t, x],x] +
(1 + Tanh[1 - b[t, x]] ) - (1 + b[t, x])*a[t, x],

(*specify the Dirichlet boundary conditions*)
bcs =
{DirichletCondition[n[t, x] == E^-t, x == 0],
DirichletCondition[b[t, x] == 1, x == 0]};

(*specify the Neumann boundary conditions (for function a(t,x))*)
BNeuma = {NeumannValue[-a[t, x], x == 0], NeumannValue[0, x == 1]};
BNeumn = NeumannValue[0, x == 1];
BNeumb = NeumannValue[0, x == 1];

(*function to use in the picewise intial conditions*)
xTild = 0.05;
n0[x_] := 1/xTild^3*(x - xTild)*(2*x^2 - xTild*x - xTild^2);
b0[x_] := 1/xTild^3 (x - xTild)*(2*x^2 - xTild*x - xTild^2);

(*initial conditions themselfs*)
incs =
{n[0, x] == Piecewise[{{n0[x], x < xTild}, {0, x > xTild}}],
a[0, x] == 0,
b[0, x] == Piecewise[{{b0[x], x < xTild}, {0, x > xTild}}]};

(*solve the system*)
{nfun, afun, bfun} =
NDSolveValue[
{op == {BNeumn, BNeuma, BNeumn}, bcs, incs},
{n, a, b}, {x, 0, 1}, {t, 0, 10}];


I get the following error:

Set::shape: Lists {nfun, afun, bfun} and NDSolveValue[<<1>>] are not the same shape.

Please, give me a hint, why my code is incorrect. If it is nessesary, I can upload a full code with all parameters, but it gives exactly the same error. My Mathematica version is 12.1.

• Looks like NDSolveValue did not return a value....(Sorry, can't check the code from a phone.) Sep 22, 2020 at 21:32
• We need to add boundary conditions at x=1. Sep 22, 2020 at 21:58
• @MichaelE2 Good point, thanks. Are any suggestions why this happens? Sep 22, 2020 at 22:17
• @AlexTrounev Thanks, I will try. But according to the tutorial reference.wolfram.com/language/FEMDocumentation/tutorial/… if no boundary conditions are specified,NDSolveValue suggestes that they are zero Neumann (which is indeed what I want according to the equations 14 from my pictures). Or it is wrong? Sep 22, 2020 at 22:19
• Normally, I would expect there to be at least one error before the Set::shape error that would given me a hint. Incorrect syntax or failure to be able to compute the solution are common reasons. Sep 22, 2020 at 23:32

In this case it is better to use not FEM orientated code as follows

eq = {-D[n[t, x], t] + D[n[t, x], x, x] -
D[n[t, x]*D[a[t, x], x], x] + a[t, x]*b[t, x] -
n[t, x] - (n[t, x])^2 ==
0, -D[a[t, x], t] +
D[a[t, x], x, x] + (1 + Tanh[1 - b[t, x]]) - (1 + b[t, x])*
a[t, x] ==
0, -D[b[t, x], t] + D[n[t, x]*D[b[t, x], x], x] - D[n[t, x], x] +
n[t, x]*D[a[t, x], x] == 0};

(*specify the Dirichlet boundary conditions*)
bcs = {n[t, x] == E^-t /. x -> 0,
b[t, x] == 1 /.
x -> 0, -Derivative[0, 1][a][t, x] + a[t, x] == 0 /. x -> 0,
Derivative[0, 1][n][t, 1] == 0, Derivative[0, 1][a][t, 1] == 0,
Derivative[0, 1][b][t, 1] == 0};

(*specify the Neumann boundary conditions (for function a(t,x))*)
BNeuma = NeumannValue[-a[t, x], x == 0];

(*function to use in the picewise intial conditions*)
xTild = 0.05;
n0[x_] := 1/xTild^3*(x - xTild)*(2*x^2 - xTild*x - xTild^2);
b0[x_] := 1/xTild^3 (x - xTild)*(2*x^2 - xTild*x - xTild^2);

(*initial conditions themselfs*)
incs = {n[0, x] == If[x <= xTild, n0[x], 0], a[0, x] == 0,
b[0, x] == If[x < xTild, b0[x], 0]};

{nfun, afun, bfun} =
NDSolveValue[Join[eq, bcs, incs], {n, a, b}, {x, 0, 1}, {t, 0, 10}];


There is message NDSolveValue::mxsst: Using maximum number of grid points 10000 allowed by the MaxPoints or MinStepSize options for independent variable x., but numerical solution exists and we can visualize it as

{DensityPlot[nfun[t, x], {x, 0, 1}, {t, 0, 10},
ColorFunction -> "Rainbow", PlotLegends -> Automatic,
FrameLabel -> Automatic, PlotRange -> All],
DensityPlot[afun[t, x], {x, 0, 1}, {t, 0, 10},
ColorFunction -> "Rainbow", PlotLegends -> Automatic,
FrameLabel -> Automatic, PlotRange -> All],
DensityPlot[bfun[t, x], {x, 0, 1}, {t, 0, 10},
ColorFunction -> "Rainbow", PlotLegends -> Automatic,
FrameLabel -> Automatic, PlotRange -> All]}


• Thank you very much! But how did you conclud that it is better to use not FEM? Is it clear from the system itself? Sep 23, 2020 at 11:45
• @АндрейКокорев It is too complicated to use FEM for system of nonlinear PDE. We know that for parabolic systems there is method of lines, So we can solve it without FEM. Sep 23, 2020 at 12:25
• thanks again, you helped a lot! Sep 23, 2020 at 13:53
• @АндрейКокорев You are welcome! Sep 23, 2020 at 13:55

As was already said, NDSolve does not evaluate. The reason is that the syntax is wrong. You should specify something like:

NDsolve[{equations, boundary conditions},..]

But you specified:

NDsolve[{{expression1,Expression2}=={ boundary conditions},..]

You must give equations, not expressions.