# Can NDSolve be used for a system with number (+40) of equations?

I have 42 coupled differential equations (6 of them are PDE's). My system has 6 objects and an object have 7 coupled equations. Since these objects are interacting with each other, their coupled equations also coupled with each other.

I can solve the system for 14 equations (2 objects). However, I can not run the code for 6 objects, NDSolve doesn't even give an error or warning. This is my first trial, and I am trying to make the simplest possible agent based model (i guess it can be referred as agent based modelling).

Is it possible to solve such a system with MMA? If so how can it be done? My final question, is there a different terminology for systems with high number of equations?

Any help would be appreciated.

The code:

Clear["Global*"]
(*SET UP*)
tmax = 2500;
(*INITIAL CONDITIONS*)

inical = SetPrecision[ 0.2001651851715455515. , 5];
iniope = 1/(1 + inical^2);
\[Mu]\[Mu] = SetPrecision[0.5500000, 5]
\[Mu]0 = SetPrecision[0.5155877999949406, 5]
(*PARAMETERS*)
Dc = 20;
kflux = 8.1;
l = 20;(*--------*)
b = 0.111 // Rationalize;
k1 = 0.7;
\[Gamma] = 2.0;
k\[Gamma] = 0.1;
\[Beta] =(*0 *)0.02;
\[Tau]n = 2.0;
k2 = 0.7;
K1 = ((kflux \[Tau]n)/k1);
\[CapitalGamma]\[CapitalGamma] = ((\[Tau]n \[Gamma] 0.9)/k1);
\[CapitalGamma] = ((\[Tau]n \[Gamma] )/k1);
T = (\[Gamma]*\[Tau]n)/k1;
K = (k\[Gamma]/k1);
K2 = k2/k1;
B = \[Beta]/k1;
D0 = Dc \[Tau]n/l^2;
Dp0 = Dp \[Tau]n /l^2;
Dp = 300;
\[Mu]1 = 0.433;
kp = 0.02;
(*q=-3.4;*)
xxx = -0.1;
qqq = -0.4;
Q1 = 0.09;
Q2 = 0.5;
ttt = 0.6;
zzz = -1.2;
sol = NDSolve[{
D[v1o1[t, x], t] ==
D0 D[v1o1[t, x], x,
x] +((\[Mu]0) K1 ((b + v1o1[t, x])/(1 + v1o1[t, x])) v2o1[t,
x]) - ((\[CapitalGamma]\[CapitalGamma] v1o1[t, x])/(K +
v1o1[t, x])) +
B -10^zzz(*Q110^xxx/\[Tau]n*)(v7o1[t, x]) v1o1[t, x],
D[v2o1[t, x], t] == 1/(1 + v1o1[t, x]^2) - v2o1[t, x],
D[v3o1[t, x],
t] == (0.001 +10^ttt (v1o1[t, x]^2/(1/k1^2 + v1o1[t, x]^2)) -
0.005*(v3o1[t, x]) -
0.001*(v3o1[t, x])*((v4o2[t, x] + v4o3[t, x])/2) +
0.1*(v5o2[t, x])),
D[v4o2[t, x],
t] == (-0.001*(v3o1[t, x])*(v4o2[t, x]) + 0.1*(v5o2[t, x]) -
0.005*(v4o2[t, x]) + 0.007),
D[v5o2[t, x], t] ==
0.001*(v3o1[t, x])*(v4o2[t, x]) - 0.1*(v5o2[t, x]) -
0.005*(v5o2[t, x]),
D[v6o2[t, x], t] == ( 0.1*(v5o2[t, x]) - 0.005*(v6o2[t, x])),
D[v7o2[t, x], t] == (
0.001 + ((0.1*(v6o2[t, x])^2)/(1 + (v6o2[t, x])^2)) -
0.005*(v7o2[t, x])),
D[v1o2[t, x], t] ==
D0 D[v1o2[t, x], x,
x] + ((\[Mu]\[Mu]) K1 ((b + v1o2[t, x])/(1 + v1o2[t, x])) v2o2[
t, x]) - ((\[CapitalGamma]\[CapitalGamma] v1o2[t, x])/(K +
v1o2[t, x])) + B -
10^zzz(*Q110^xxx/\[Tau]n*)(v7o2[t, x]) v1o2[t, x],
D[v2o2[t, x], t] == 1/(1 + v1o2[t, x]^2) - v2o2[t, x],
D[v3o2[t, x],
t] == (0.001 +10^ttt (v1o2[t, x]^2/(1/k1^2 + v1o2[t, x]^2)) -
0.005*(v3o2[t, x]) -
0.001*(v3o2[t, x])*((v4o1[t, x] + v4o3[t, x])/2) +
0.1*(v5o1[t, x])),
D[v4o1[t, x],
t] == (-0.001*(v3o2[t, x])*(v4o1[t, x]) + 0.1*(v5o1[t, x]) -
0.005*(v4o1[t, x]) + 0.007),
D[v5o1[t, x], t] ==
0.001*(v3o2[t, x])*(v4o1[t, x]) - 0.1*(v5o1[t, x]) -
0.005*(v5o1[t, x]),
D[v6o1[t, x], t] == ( 0.1*(v5o1[t, x]) - 0.005*(v6o1[t, x])),
D[v7o1[t, x], t] == (
0.001 + ((0.1*(v6o1[t, x])^2)/(1 + (v6o1[t, x])^2)) -
0.005*(v7o1[t, x])),
(**)
D[v1o3[t, x], t] ==
D0 D[v1o3[t, x], x,
x] + ((\[Mu]0) K1 ((b + v1o3[t, x])/(1 + v1o3[t, x])) v2o3[t,
x]) - ((\[CapitalGamma]\[CapitalGamma] v1o3[t, x])/(K +
v1o3[t, x])) + B -
10^zzz(*Q110^xxx/\[Tau]n*)(v7o3[t, x]) v1o3[t, x],
D[v2o3[t, x], t] == 1/(1 + v1o3[t, x]^2) - v2o3[t, x],
D[v3o3[t, x],
t] == (0.001 +10^ttt (v1o3[t, x]^2/(1/k1^2 + v1o3[t, x]^2)) -
0.005*(v3o3[t, x]) -
0.001*(v3o3[t, x])*((v4o1[t, x] + v4o4[t, x])/2) +
0.1*(v5o4[t, x])),
D[v4o4[t, x],
t] == (-0.001*(v3o3[t, x])*(v4o4[t, x]) + 0.1*(v5o4[t, x]) -
0.005*(v4o4[t, x]) + 0.007),
D[v5o4[t, x], t] ==
0.001*(v3o3[t, x])*(v4o4[t, x]) - 0.1*(v5o4[t, x]) -
0.005*(v5o4[t, x]),
D[v6o4[t, x], t] == ( 0.1*(v5o4[t, x]) - 0.005*(v6o4[t, x])),
D[v7o4[t, x], t] == (
0.001 + ((0.1*(v6o4[t, x])^2)/(1 + (v6o4[t, x])^2)) -
0.005*(v7o4[t, x])),
D[v1o4[t, x], t] ==
D0 D[v1o4[t, x], x,
x] + ((\[Mu]\[Mu]) K1 ((b + v1o4[t, x])/(1 + v1o4[t, x])) v2o4[
t, x]) - ((\[CapitalGamma]\[CapitalGamma] v1o4[t, x])/(K +
v1o4[t, x])) + B -
10^zzz(*Q110^xxx/\[Tau]n*)(v7o4[t, x]) v1o4[t, x],
D[v2o4[t, x], t] == 1/(1 + v1o4[t, x]^2) - v2o4[t, x],
D[v3o4[t, x],
t] == (0.001 +10^ttt (v1o4[t, x]^2/(1/k1^2 + v1o4[t, x]^2)) -
0.005*(v3o4[t, x]) -
0.001*(v3o4[t, x])*((v4o1[t, x] + v4o3[t, x])/2) +
0.1*(v5o3[t, x])),
D[v4o3[t, x],
t] == (-0.001*(v3o4[t, x])*(v4o3[t, x]) + 0.1*(v5o3[t, x]) -
0.005*(v4o3[t, x]) + 0.007),
D[v5o3[t, x], t] ==
0.001*(v3o4[t, x])*(v4o3[t, x]) - 0.1*(v5o3[t, x]) -
0.005*(v5o3[t, x]),
D[v6o3[t, x], t] == ( 0.1*(v5o3[t, x]) - 0.005*(v6o3[t, x])),
D[v7o3[t, x], t] == (
0.001 + ((0.1*(v6o3[t, x])^2)/(1 + (v6o3[t, x])^2)) -
0.005*(v7o3[t, x])),
(**)
D[v1o5[t, x], t] ==
D0 D[v1o5[t, x], x,
x] + ((\[Mu]0) K1 ((b + v1o5[t, x])/(1 + v1o5[t, x])) v2o5[t,
x]) - ((\[CapitalGamma]\[CapitalGamma] v1o5[t, x])/(K +
v1o5[t, x])) + B -
10^zzz(*Q110^xxx/\[Tau]n*)(v7o5[t, x]) v1o5[t, x],
D[v2o5[t, x], t] == 1/(1 + v1o5[t, x]^2) - v2o5[t, x],
D[v3o5[t, x],
t] == (0.001 +10^ttt (v1o5[t, x]^2/(1/k1^2 + v1o5[t, x]^2)) -
0.005*(v3o5[t, x]) -
0.001*(v3o5[t, x])*((v4o6[t, x] + v4o3[t, x])/2) +
0.1*(v5o6[t, x])),
D[v4o6[t, x],
t] == (-0.001*(v3o5[t, x])*(v4o6[t, x]) + 0.1*(v5o6[t, x]) -
0.005*(v4o6[t, x]) + 0.007),
D[v5o6[t, x], t] ==
0.001*(v3o5[t, x])*(v4o6[t, x]) - 0.1*(v5o6[t, x]) -
0.005*(v5o6[t, x]),
D[v6 o6[t, x], t] == ( 0.1*(v5o6[t, x]) - 0.005*(v6o6[t, x])),
D[v7o6[t, x], t] == (
0.001 + ((0.1*(v6o6[t, x])^2)/(1 + (v6o6[t, x])^2)) -
0.005*(v7o6[t, x])),
D[v1o6[t, x], t] ==
D0 D[v1o6[t, x], x,
x] + ((\[Mu]\[Mu]) K1 ((b + v1o6[t, x])/(1 + v1o6[t, x])) v2o6[
t, x]) - ((\[CapitalGamma]\[CapitalGamma] v1o6[t, x])/(K +
v1o6[t, x])) + B -
10^zzz(*Q110^xxx/\[Tau]n*)(v7o6[t, x]) v1o6[t, x],
D[v2o6[t, x], t] == 1/(1 + v1o6[t, x]^2) - v2o6[t, x],
D[v3o6[t, x],
t] == (0.001 +10^ttt (v1o6[t, x]^2/(1/k1^2 + v1o6[t, x]^2)) -
0.005*(v3o6[t, x]) -
0.001*(v3o6[t, x])*((v4o5[t, x] + v4o3[t, x])/2) +
0.1*(v5o5[t, x])),
D[v4o5[t, x],
t] == (-0.001*(v3o6[t, x])*(v4o5[t, x]) + 0.1*(v5o5[t, x]) -
0.005*(v4o5[t, x]) + 0.007),
D[v5o5[t, x], t] ==
0.001*(v3o6[t, x])*(v4o5[t, x]) - 0.1*(v5o5[t, x]) -
0.005*(v5o5[t, x]),
D[v6o5[t, x], t] == ( 0.1*(v5o5[t, x]) - 0.005*(v6o5[t, x])),
D[v7o5[t, x], t] == (
0.001 + ((0.1*(v6o5[t, x])^2)/(1 + (v6o5[t, x])^2)) -
0.005*(v7o5[t, x])),
v2o1[0, x] == iniope,
v1o1[0, x] == If[-20 < x < 20, 2 , inical],
v1o1[t, -200] == inical ,
v1o1[t, 200] == inical ,
v3o1[0, x] == 0,
v4o1[0, x] == 0.1,(**)
v5o1[0, x] == 0,
v6o1[0, x] == 0,
v7o1[0, x] == 0,
v2o2[0, x] == iniope,
v1o2[0, x] == If[-20 < x < 20, 2 , inical],
v1o2[t, -200] == inical ,
v1o2[t, 200] == inical ,
v3o2[0, x] == 0,
v4o2[0, x] == 0.1,(**)
v5o2[0, x] == 0,
v6o2[0, x] == 0,
v7o2[0, x] == 0,
(**)
v2o3[0, x] == iniope,
v1o3[0, x] == If[-20 < x < 20, 2 , inical],
v1o3[t, -200] == inical ,
v1o3[t, 200] == inical ,
v3o3[0, x] == 0,
v4o3[0, x] == 0.1,(**)
v5o3[0, x] == 0,
v6o3[0, x] == 0,
v7o3[0, x] == 0,
v2o4[0, x] == iniope,
v1o4[0, x] == If[-20 < x < 20, 2 , inical],
v1o4[t, -200] == inical ,
v1o4[t, 200] == inical ,
v3o4[0, x] == 0,
v4o4[0, x] == 0.1,(**)
v5o4[0, x] == 0,
v6o4[0, x] == 0,
v7o4[0, x] == 0,
(**)
v2o5[0, x] == iniope,
v1o5[0, x] == If[-20 < x < 20, 2 , inical],
v1o5[t, -200] == inical ,
v1o5[t, 200] == inical ,
v3o5[0, x] == 0,
v4o5[0, x] == 0.1,(**)
v5o5[0, x] == 0,
v6o5[0, x] == 0,
v7o5[0, x] == 0,
v2o6[0, x] == iniope,
v1o6[0, x] == If[-20 < x < 20, 2 , inical],
v1o6[t, -200] == inical ,
v1o6[t, 200] == inical ,
v3o6[0, x] == 0,
v4o6[0, x] == 0.1,(**)
v5o6[0, x] == 0,
v6o6[0, x] == 0,
v7o6[0, x] == 0
},
{v1o1, v2o1, v3o1, v4o1, v5o1, v6o1, v7o1, v1o2, v2o2, v3o2, v4o2,
v5o2, v6o2, v7o2(**), v1o3, v2o3, v3o3, v4o3, v5o3, v6o3, v7o3,
v1o4, v2o4, v3o4, v4o4, v5o4, v6o4, v7o4,(**)v1o5, v2o5, v3o5, v4o5,
v5o5, v6o5, v7o5, v1o6, v2o6, v3o6, v4o6, v5o6, v6o6, v7o6}, {t, 0,
tmax}, {x, -200, 200},
Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid",
"MaxPoints" -> 101}}]
Plot[{Evaluate[{v3o1[t, 0]} /. sol], Evaluate[{v3o2[t, 0]} /. sol],
Evaluate[{v3o3[t, 0]} /. sol], Evaluate[{v3o4[t, 0]} /. sol],
Evaluate[{v3o5[t, 0]} /. sol], Evaluate[{v3o6[t, 0]} /. sol]}, {t,
0, tmax}, PlotStyle -> {Red, Blue}, ColorFunction -> "Rainbow",
PlotLegends -> Automatic]
$$$$

• in theory, with a numerical solver, having 16 equations vs. 32 vs. 64, should not make any difference (except for memory and time). As to why NDSolve can't solve your system when you add more equations, it is hard to know, since you did not show an example or say what error or what happens exactly. does it hang? what does it do? Without an example, hard to say. Commented Feb 28, 2021 at 9:28
• @Nasser you are right, I added the code. Commented Feb 28, 2021 at 10:15
• @AlexTrounev That's a PDE, isn't it? An ODE has the form u'[t] == F[t, u[t]....I think we should be strict in our definitions when dealing with a computer program, or at least as strict as the program. Commented Feb 28, 2021 at 15:00
• @confused Quitting the kernel is known to fix inexplicable problems. If it works, it should be because of a forgotten definition that wasn't cleared. Commented Feb 28, 2021 at 15:01
• Quitting the kernel could also reset a bug. Like rebooting the computer. :) Commented Feb 28, 2021 at 15:15

After some debugging this code works on v12.0 and 12.2 as well

Clear["Global*"]
(*SET UP*)
tmax = 2500;
(*INITIAL CONDITIONS*)

inical = SetPrecision[0.2001651851715455515., 5];
iniope = 1/(1 + inical^2);
\[Mu]\[Mu] = SetPrecision[0.5500000, 5];
\[Mu]0 = SetPrecision[0.5155877999949406, 5];
(*PARAMETERS*)
Dc = 20;
kflux = 8.1;
l = 20;(*--------*)b = 0.111 // Rationalize;
k1 = 0.7;
\[Gamma] = 2.0;
k\[Gamma] = 0.1;
\[Beta] =(*0*)0.02;
\[Tau]n = 2.0;
k2 = 0.7;
K1 = ((kflux \[Tau]n)/k1);
\[CapitalGamma]\[CapitalGamma] = ((\[Tau]n \[Gamma] 0.9)/k1);
\[CapitalGamma] = ((\[Tau]n \[Gamma])/k1);
T = (\[Gamma]*\[Tau]n)/k1;
K = (k\[Gamma]/k1);
K2 = k2/k1;
B = \[Beta]/k1;
D0 = Dc \[Tau]n/l^2;
Dp0 = Dp \[Tau]n/l^2;
Dp = 300;
\[Mu]1 = 0.433;
kp = 0.02;
(*q=-3.4;*)
xxx = -0.1;
qqq = -0.4;
Q1 = 0.09;
Q2 = 0.5;
ttt = 0.6;
zzz = -1.2; var = {v1o1, v2o1, v3o1, v4o1, v5o1, v6o1, v7o1, v1o2,
v2o2, v3o2, v4o2, v5o2, v6o2, v7o2, v1o3, v2o3, v3o3, v4o3, v5o3,
v6o3, v7o3, v1o4, v2o4, v3o4, v4o4, v5o4, v6o4, v7o4, v1o5, v2o5,
v3o5, v4o5, v5o5, v6o5, v7o5, v1o6, v2o6, v3o6, v4o6, v5o6, v6o6,
v7o6};
eq = {D[v1o1[t, x], t] ==
D0 D[v1o1[t, x], x,
x] + ((\[Mu]0) K1 ((b + v1o1[t, x])/(1 + v1o1[t, x])) v2o1[t,
x]) - ((\[CapitalGamma]\[CapitalGamma] v1o1[t, x])/(K +
v1o1[t, x])) + B - 10^zzz (v7o1[t, x]) v1o1[t, x],
D[v2o1[t, x], t] == 1/(1 + v1o1[t, x]^2) - v2o1[t, x],
D[v3o1[t, x],
t] == (0.001 + 10^ttt (v1o1[t, x]^2/(1/k1^2 + v1o1[t, x]^2)) -
0.005*(v3o1[t, x]) -
0.001*(v3o1[t, x])*((v4o2[t, x] + v4o3[t, x])/2) +
0.1*(v5o2[t, x])),
D[v4o2[t, x],
t] == (-0.001*(v3o1[t, x])*(v4o2[t, x]) + 0.1*(v5o2[t, x]) -
0.005*(v4o2[t, x]) + 0.007),
D[v5o2[t, x], t] ==
0.001*(v3o1[t, x])*(v4o2[t, x]) - 0.1*(v5o2[t, x]) -
0.005*(v5o2[t, x]),
D[v6o2[t, x], t] == (0.1*(v5o2[t, x]) - 0.005*(v6o2[t, x])),
D[v7o2[t, x],
t] == (0.001 + ((0.1*(v6o2[t, x])^2)/(1 + (v6o2[t, x])^2)) -
0.005*(v7o2[t, x])),
D[v1o2[t, x], t] ==
D0 D[v1o2[t, x], x,
x] + ((\[Mu]\[Mu]) K1 ((b + v1o2[t, x])/(1 + v1o2[t, x])) v2o2[
t, x]) - ((\[CapitalGamma]\[CapitalGamma] v1o2[t, x])/(K +
v1o2[t, x])) + B - 10^zzz (v7o2[t, x]) v1o2[t, x],
D[v2o2[t, x], t] == 1/(1 + v1o2[t, x]^2) - v2o2[t, x],
D[v3o2[t, x],
t] == (0.001 + 10^ttt (v1o2[t, x]^2/(1/k1^2 + v1o2[t, x]^2)) -
0.005*(v3o2[t, x]) -
0.001*(v3o2[t, x])*((v4o1[t, x] + v4o3[t, x])/2) +
0.1*(v5o1[t, x])),
D[v4o1[t, x],
t] == (-0.001*(v3o2[t, x])*(v4o1[t, x]) + 0.1*(v5o1[t, x]) -
0.005*(v4o1[t, x]) + 0.007),
D[v5o1[t, x], t] ==
0.001*(v3o2[t, x])*(v4o1[t, x]) - 0.1*(v5o1[t, x]) -
0.005*(v5o1[t, x]),
D[v6o1[t, x], t] == (0.1*(v5o1[t, x]) - 0.005*(v6o1[t, x])),
D[v7o1[t, x],
t] == (0.001 + ((0.1*(v6o1[t, x])^2)/(1 + (v6o1[t, x])^2)) -
0.005*(v7o1[t, x])),(**)
D[v1o3[t, x], t] ==
D0 D[v1o3[t, x], x,
x] + ((\[Mu]0) K1 ((b + v1o3[t, x])/(1 + v1o3[t, x])) v2o3[t,
x]) - ((\[CapitalGamma]\[CapitalGamma] v1o3[t, x])/(K +
v1o3[t, x])) + B - 10^zzz (v7o3[t, x]) v1o3[t, x],
D[v2o3[t, x], t] == 1/(1 + v1o3[t, x]^2) - v2o3[t, x],
D[v3o3[t, x],
t] == (0.001 + 10^ttt (v1o3[t, x]^2/(1/k1^2 + v1o3[t, x]^2)) -
0.005*(v3o3[t, x]) -
0.001*(v3o3[t, x])*((v4o1[t, x] + v4o4[t, x])/2) +
0.1*(v5o4[t, x])),
D[v4o4[t, x],
t] == (-0.001*(v3o3[t, x])*(v4o4[t, x]) + 0.1*(v5o4[t, x]) -
0.005*(v4o4[t, x]) + 0.007),
D[v5o4[t, x], t] ==
0.001*(v3o3[t, x])*(v4o4[t, x]) - 0.1*(v5o4[t, x]) -
0.005*(v5o4[t, x]),
D[v6o4[t, x], t] == (0.1*(v5o4[t, x]) - 0.005*(v6o4[t, x])),
D[v7o4[t, x],
t] == (0.001 + ((0.1*(v6o4[t, x])^2)/(1 + (v6o4[t, x])^2)) -
0.005*(v7o4[t, x])),
D[v1o4[t, x], t] ==
D0 D[v1o4[t, x], x,
x] + ((\[Mu]\[Mu]) K1 ((b + v1o4[t, x])/(1 + v1o4[t, x])) v2o4[
t, x]) - ((\[CapitalGamma]\[CapitalGamma] v1o4[t, x])/(K +
v1o4[t, x])) + B - 10^zzz (v7o4[t, x]) v1o4[t, x],
D[v2o4[t, x], t] == 1/(1 + v1o4[t, x]^2) - v2o4[t, x],
D[v3o4[t, x],
t] == (0.001 + 10^ttt (v1o4[t, x]^2/(1/k1^2 + v1o4[t, x]^2)) -
0.005*(v3o4[t, x]) -
0.001*(v3o4[t, x])*((v4o1[t, x] + v4o3[t, x])/2) +
0.1*(v5o3[t, x])),
D[v4o3[t, x],
t] == (-0.001*(v3o4[t, x])*(v4o3[t, x]) + 0.1*(v5o3[t, x]) -
0.005*(v4o3[t, x]) + 0.007),
D[v5o3[t, x], t] ==
0.001*(v3o4[t, x])*(v4o3[t, x]) - 0.1*(v5o3[t, x]) -
0.005*(v5o3[t, x]),
D[v6o3[t, x], t] == (0.1*(v5o3[t, x]) - 0.005*(v6o3[t, x])),
D[v7o3[t, x],
t] == (0.001 + ((0.1*(v6o3[t, x])^2)/(1 + (v6o3[t, x])^2)) -
0.005*(v7o3[t, x])),(**)
D[v1o5[t, x], t] ==
D0 D[v1o5[t, x], x,
x] + ((\[Mu]0) K1 ((b + v1o5[t, x])/(1 + v1o5[t, x])) v2o5[t,
x]) - ((\[CapitalGamma]\[CapitalGamma] v1o5[t, x])/(K +
v1o5[t, x])) + B - 10^zzz (v7o5[t, x]) v1o5[t, x],
D[v2o5[t, x], t] == 1/(1 + v1o5[t, x]^2) - v2o5[t, x],
D[v3o5[t, x],
t] == (0.001 + 10^ttt (v1o5[t, x]^2/(1/k1^2 + v1o5[t, x]^2)) -
0.005*(v3o5[t, x]) -
0.001*(v3o5[t, x])*((v4o6[t, x] + v4o3[t, x])/2) +
0.1*(v5o6[t, x])),
D[v4o6[t, x],
t] == (-0.001*(v3o5[t, x])*(v4o6[t, x]) + 0.1*(v5o6[t, x]) -
0.005*(v4o6[t, x]) + 0.007),
D[v5o6[t, x], t] ==
0.001*(v3o5[t, x])*(v4o6[t, x]) - 0.1*(v5o6[t, x]) -
0.005*(v5o6[t, x]),
D[v6o6[t, x], t] == (0.1*(v5o6[t, x]) - 0.005*(v6o6[t, x])),
D[v7o6[t, x],
t] == (0.001 + ((0.1*(v6o6[t, x])^2)/(1 + (v6o6[t, x])^2)) -
0.005*(v7o6[t, x])),
D[v1o6[t, x], t] ==
D0 D[v1o6[t, x], x,
x] + ((\[Mu]\[Mu]) K1 ((b + v1o6[t, x])/(1 + v1o6[t, x])) v2o6[
t, x]) - ((\[CapitalGamma]\[CapitalGamma] v1o6[t, x])/(K +
v1o6[t, x])) + B - 10^zzz (v7o6[t, x]) v1o6[t, x],
D[v2o6[t, x], t] == 1/(1 + v1o6[t, x]^2) - v2o6[t, x],
D[v3o6[t, x],
t] == (0.001 + 10^ttt (v1o6[t, x]^2/(1/k1^2 + v1o6[t, x]^2)) -
0.005*(v3o6[t, x]) -
0.001*(v3o6[t, x])*((v4o5[t, x] + v4o3[t, x])/2) +
0.1*(v5o5[t, x])),
D[v4o5[t, x],
t] == (-0.001*(v3o6[t, x])*(v4o5[t, x]) + 0.1*(v5o5[t, x]) -
0.005*(v4o5[t, x]) + 0.007),
D[v5o5[t, x], t] ==
0.001*(v3o6[t, x])*(v4o5[t, x]) - 0.1*(v5o5[t, x]) -
0.005*(v5o5[t, x]),
D[v6o5[t, x], t] == (0.1*(v5o5[t, x]) - 0.005*(v6o5[t, x])),
D[v7o5[t, x],
t] == (0.001 + ((0.1*(v6o5[t, x])^2)/(1 + (v6o5[t, x])^2)) -
0.005*(v7o5[t, x]))};
icbc = {v2o1[0, x] == iniope,
v1o1[0, x] == If[-20 < x < 20, 2, inical], v1o1[t, -200] == inical,
v1o1[t, 200] == inical, v3o1[0, x] == 0, v4o1[0, x] == 0.1,(**)
v5o1[0, x] == 0, v6o1[0, x] == 0, v7o1[0, x] == 0,
v2o2[0, x] == iniope, v1o2[0, x] == If[-20 < x < 20, 2, inical],
v1o2[t, -200] == inical, v1o2[t, 200] == inical, v3o2[0, x] == 0,
v4o2[0, x] == 0.1,(**)v5o2[0, x] == 0, v6o2[0, x] == 0,
v7o2[0, x] == 0,(**)v2o3[0, x] == iniope,
v1o3[0, x] == If[-20 < x < 20, 2, inical], v1o3[t, -200] == inical,
v1o3[t, 200] == inical, v3o3[0, x] == 0, v4o3[0, x] == 0.1,(**)
v5o3[0, x] == 0, v6o3[0, x] == 0, v7o3[0, x] == 0,
v2o4[0, x] == iniope, v1o4[0, x] == If[-20 < x < 20, 2, inical],
v1o4[t, -200] == inical, v1o4[t, 200] == inical, v3o4[0, x] == 0,
v4o4[0, x] == 0.1,(**)v5o4[0, x] == 0, v6o4[0, x] == 0,
v7o4[0, x] == 0,(**)v2o5[0, x] == iniope,
v1o5[0, x] == If[-20 < x < 20, 2, inical], v1o5[t, -200] == inical,
v1o5[t, 200] == inical, v3o5[0, x] == 0, v4o5[0, x] == 0.1,(**)
v5o5[0, x] == 0, v6o5[0, x] == 0, v7o5[0, x] == 0,
v2o6[0, x] == iniope, v1o6[0, x] == If[-20 < x < 20, 2, inical],
v1o6[t, -200] == inical, v1o6[t, 200] == inical, v3o6[0, x] == 0,
v4o6[0, x] == 0.1,(**)v5o6[0, x] == 0, v6o6[0, x] == 0,
v7o6[0, x] == 0};
sol = NDSolve[{eq, icbc}, var, {t, 0, tmax}, {x, -200, 200},
Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid",
"MaxPoints" -> 101}}];


Visualization of numerical solution

Table[DensityPlot[
Evaluate[var[[i]][t, x] /. sol[[1]]], {t, 0, tmax}, {x, -200, 200},
PlotLabel -> var[[i]], ColorFunction -> "Rainbow",
PlotLegends -> Automatic, PlotRange -> All], {i, Length[var]}]
`