# Trying to use NDSolve to solve a system of non-linear PDEs representing the development of Alzheimer's disease over time

I'm trying to do something very similar to this question, except with a much larger system of PDEs from this paper. The PDEs are equations 1-18 in that paper. Equations 9, 10 and 15-18 have laplacian terms in them. From the paper, I know that a rectangular domain is used and that the variables that have laplacian terms in their equations satisfy periodic boundary conditions. Therefore, I wrote the following code:

    (* Parameters *)
DAO = 4.32*10^-2;
DH = 8.11*10^-2;
DTAlpha = 6.55*10^-2;
DTBeta = 6.55*10^-2;
Di10 = 6.04*10^-2;
DP = 1.2*10^-1;
lambdaBetaI = 9.51*10^-6;
lambdaN = 8*10^-9;
lambdaA = 8*10^-10;
lambdaTau0 = 8.1*10^-11;
lambdaTau = 1.35*10^-11;
lambdaF = 1.662*10^-3;
lambdaATAlpha = 1.54;
lambdaAABeta0 = 1.793;
lambdaAO = 5*10^-2;
lambdaH = 3*10^-5;
lambdaMF = 2*10^-2;
lambdaMA = 2.3*10^-3;
lambdaM1TBeta = 6*10^-3;
lambdaM1HatTBeta = 6*10^-4;
lambdaTBetaM2 = 1.5 * 10^-2;
lambdaTBetaM2Hat = 1.5 * 10^-2;
lambdaTAlphaM1 = 3 * 10^-2;
lambdaTAlphaM1Hat = 3 * 10^-2;
lambdai10M2 = 6.67 * 10^-3;
lambdai10M2Hat = 6.67 * 10^-3;
lambdaPA = 6.6 * 10^-8;
lambdaPM2 = 1.32 * 10^-7;
myTheta = 0.9;
alpha1 = 5;
myBeta = 10;
myGamma = 1;
dABetaI = 9.51;
dABeta0M = 9.51;
dABeta0M = 2 * 10^-3;
dABeta0MHat = 10^-2;
dTau = 0.277;
dFI = 2.77 * 10^-3;
dFO = 2.77 * 10^-4;
dN = 1.9 * 10^-4;
dNF = 3.4 * 10^-4;
dNT = 1.7 * 10^-4;
dNdM = 0.06;
dNdMhat = 0.02;
dA = 1.2 * 10^-3;
dM1 = 0.015;
dM2 = 0.015;
dM1Hat = 0.015;
dM2Hat = 0.015;
dAO = 0.951;
dH = 58.71;
dTAlpha = 55.45;
dTBeta = 3.33 * 10^2;
di10 = 16.64;
dP = 1.73;
r0 = 6;
m0 = 5 * 10^-2;
n0 = 0.14;
mG0 = 0.047;
a0 = 0.14;
kBarABeta0 = 7 * 10^-3;
kBarNd = 10^-3;
ki10 = 2.5 * 10^-6;
kTBeta = 2.5 * 10^-7;
kM = 0.047;
kMHat = 0.047;
kM1 = 0.03;
kM2 = 0.017;
kM1Hat = 0.04;
kM2Hat = 0.007;
kfI = 3.36 * 10^-10;
kFO = 2.58 * 10^-11;
kAO = 10^-7;
kP = 6 * 10^-9;
ktAlpha = 4 * 10^-5;

(* System of Nonlinear PDEs *)
pde ={D[aBetaI[x,y,t],t] == (lambdaBetaI* (1 +If[0<=t<=100,r0 * t / 100,r0]) - dABetaI * aBetaI[x,y,t]) * (n[x,y,t] / n0),
D[aBeta0[x,y,t],t] == aBetaI[x,y,t] * Abs[D[n[x,y,t],t]] + lambdaN * (n[x,y,t] / n0) + lambdaA * (a[x,y,t] / a0) - (dABeta0MHat * (m1Hat[x,y,t] + myTheta * m2Hat[x,y,t]) + dABeta0M * (m1[x,y,t] + myTheta * m2[x,y,t])) * (aBeta0[x,y,t] / (aBeta0[x,y,t] + kBarABeta0)),
D[myTau[x,y,t],t] == (lambdaTau0 + lambdaTau *If[0<=t<=100,r0 * t / 100,r0] - dTau * myTau[x,y,t]) * (n[x,y,t] / n0),
D[fI[x,y,t],t] == (lambdaF * myTau[x,y,t] - dFI * fI[x,y,t])*(n[x,y,t] / n0),
D[fO[x,y,t],t] == fI[x,y,t] * Abs[D[n[x,y,t],t]] - dFO * fO[x,y,t],
D[n[x,y,t],t] == -dNF * (fI[x,y,t] / (fI[x,y,t] + kfI)) * n[x,y,t] - dNT * (tAlpha[x,y,t] / (tAlpha[x,y,t] + ktAlpha)) * (1 / (1 + myGamma * (i10[x,y,t] / ki10)))*n[x,y,t],
D[a[x,y,t],t] == lambdaAABeta0 * aBeta0[x,y,t] + lambdaATAlpha * tAlpha[x,y,t] - dA * a[x,y,t],
D[nD[x,y,t],t] == dNF * (fI[x,y,t] / (fI[x,y,t] + kfI)) * n[x,y,t] + dNT * (tAlpha[x,y,t] / (tAlpha[x,y,t] + ktAlpha)) * (1 / (1 + myGamma * (i10[x,y,t] / ki10)))*n[x,y,t] - dNdM * (m1[x,y,t] + m2[x,y,t]) * (nD[x,y,t] / (nD[x,y,t] + kBarNd)) - dNdMhat * (m1Hat[x,y,t] + m2Hat[x,y,t])*(nD[x,y,t] / (nD[x,y,t] + kBarNd)),
D[aO[x,y,t],t] == lambdaAO * aBeta0[x,y,t] - dAO * aO[x,y,t] + DAO * Laplacian[aO[x,y,t],{x,y,t}],
D[h[x,y,t],t] == lambdaH *nD[x,y,t] - dH * h[x,y,t] + DH * Laplacian[h[x,y,t],{x,y,t}],
D[m1[x,y,t],t] == mG0 * (lambdaMF * (fO[x,y,t] / (fO[x,y,t] + kFO)) + lambdaMA * (aO[x,y,t] / (aO[x,y,t] + kAO))) * ((myBeta * (tAlpha[x,y,t] / (tAlpha[x,y,t] + ktAlpha))) /(myBeta *(tAlpha[x,y,t] / (tAlpha[x,y,t] + ktAlpha)) +(i10[x,y,t] / (i10[x,y,t] + ki10)))) - lambdaM1TBeta * (tBeta[x,y,t] / (tBeta[x,y,t] + kTBeta)) * m1[x,y,t] - dM1 * m1[x,y,t] - Div[m1[x,y,t] * Grad[h[x,y,t],{x,y,t}],{x,y,t}],
D[m2[x,y,t],t] == mG0 * (lambdaMF * (fO[x,y,t] / (fO[x,y,t] + kFO)) + lambdaMA * (aO[x,y,t] / (aO[x,y,t] + kAO))) * ((i10[x,y,t] / (i10[x,y,t] + kI10)) /(myBeta * (tAlpha[x,y,t] / (tAlpha[x,y,t] + ktAlpha)) + (i10[x,y,t] / (i10[x,y,t] + kI10)))) + lambdaM1TBeta * (tBeta[x,y,t] / (tBeta[x,y,t] + kTBeta)) * m1[x,y,t] - dM2 * m2[x,y,t] - Div[m2[x,y,t] * Grad[h[x,y,t],{x,y,t}],{x,y,t}],
D[m1Hat[x,y,t],t] == (alpha1 * (p[x,y,t] / (p[x,y,t] + kP))) * (m0 - ( m1Hat[x,y,t] + m2Hat[x,y,t])) *  ((myBeta * (tAlpha[x,y,t] / (tAlpha[x,y,t] + ktAlpha))) /(myBeta * (tAlpha[x,y,t] / (tAlpha[x,y,t] + ktAlpha)) + (i10[x,y,t] / (i10[x,y,t] + ki10)))) - lambdaM1HatTBeta * (tBeta[x,y,t] / (tBeta[x,y,t] + kTBeta)) * m1Hat[x,y,t] - dM1Hat * m1Hat[x,y,t] - Div[m1Hat[x,y,t] * Grad[aO[x,y,t],{x,y,t}],{x,y,t}],
D[m2Hat[x,y,t],t] == (alpha1 * (p[x,y,t] / (p[x,y,t] + kP))) * (m0 - ( m1Hat[x,y,t] + m2Hat[x,y,t])) *  ((i10[x,y,t] / (i10[x,y,t] + ki10)) /(myBeta * (tAlpha[x,y,t] / (tAlpha[x,y,t] + ktAlpha)) + (i10[x,y,t] / (i10[x,y,t] + ki10)))) + lambdaM1HatTBeta * (tBeta[x,y,t] / (tBeta[x,y,t] + kTBeta)) * m1Hat[x,y,t] - dM2Hat * m2Hat[x,y,t] - Div[m2Hat[x,y,t] * Grad[aO[x,y,t],{x,y,t}],{x,y,t}],
D[tBeta[x,y,t],t] == lambdaTBetaM2 *m2[x,y,t]  + lambdaTBetaM2Hat * m2Hat[x,y,t] - dTBeta * tBeta[x,y,t] + DTBeta * Laplacian[tBeta[x,y,t],{x,y,t}],
D[tAlpha[x,y,t],t] == lambdaTAlphaM1 *m1[x,y,t]  + lambdaTAlphaM1Hat * m1Hat[x,y,t] - dTAlpha * tAlpha[x,y,t] + DTAlpha * Laplacian[tAlpha[x,y,t],{x,y,t}],
D[i10[x,y,t],t] == lambdai10M2 *m2[x,y,t]  + lambdai10M2Hat * m2Hat[x,y,t] - di10 * i10[x,y,t] + Di10 * Laplacian[i10[x,y,t],{x,y,t}],
D[p[x,y,t],t]== lambdaPA *a[x,y,t]  + lambdaPM2 * m2[x,y,t] - dP * p[x,y,t] + DP * Laplacian[p[x,y,t],{x,y,t}]} ;

(* Periodic Boundary Conditions *)
bc = {aO[0,y,t] == aO[1,y,t], aO[x,0,t] == aO[x,1,t],h[0,y,t] ==h[1,y,t], h[x,0,t] == h[x,1,t], tBeta[0,y,t] == tBeta[1,y,t], tBeta[x,0,t] == tBeta[x,1,t], i10[0,y,t] == i10[1,y,t], i10[x,0,t] == i10[x,1,t], tAlpha[0,y,t] == tAlpha[1,y,t], tAlpha[x,0,t] == tAlpha[x,1,t], p[0,y,t] == p[1,y,t], p[x,0,t] == p[x,1,t]};

(* Initial Conditions *)
ic = {aBetaI[x,y,0] == 10^(-6),aBeta0[x,y,0] == 10^(-8),myTau[x,y,0] == 1.37*10^(-10),fI[x,y,0] == 3.36*10^(-10),fO[x,y,0] == 3.36*10^(-11),n[x,y,0] == 0.14,a[x,y,0] == 0.14,m1[x,y,0] == 0.02,m2[x,y,0] == 0.02,m1Hat[x,y,0] == 0,m2Hat[x,y,0] == 0, nD[x,y,0] == 0,h[x,y,0] == 1.3*10^(-11),tBeta[x,y,0] == 10^(-6),tAlpha[x,y,0] == 2*10^{-5},i10[x,y,0] == 10^(-5),p[x,y,0] == 5*10^(-9),aO[x,y,0] == 0};

eqns = Flatten@{pde, bc, ic};

sols = NDSolve[eqns,{aBetaI,aBeta0,myTau,fI,fO,n,a,m1,m2,m1Hat,m2Hat,nD,h,tBeta,tAlpha,i10,p,aO},{x,0,1},{y,0,1},{t,0,3650}, Method -> {"MethodOfLines", "SpatialDiscretization" -> {"TensorProductGrid", "MaxPoints" -> 200}}];


I thought that I should be able to use the same method that is used here but I get the following error:

NDSolve: Cannot solve to find an explicit formula for the derivatives. Consider using the option Method -> {"EquationSimplification"->"Residual"}.


Does anyone see any reason why I can't use the method that I'm trying to use? Also, should the variables that do not have laplacian terms only be a function of t? Currently, every variable is a function of x, y and t.

UPDATE: When I replace the last line with the following:

sols = NDSolve[eqns {aBetaI,aBeta0,myTau,fI,fO,n,a,m1,m2,m1Hat,m2Hat,nD,h,tBeta,tAlpha,i10,p,aO},{x,0,1},{y,0,1},{t,0,3650},Method ->{"EquationSimplification"->"Residual"}]


I get an error:

• Are you sure you wrote the system of equations 1-18 correctly? – Alex Trounev Oct 16 '19 at 16:13
• I looked over them a couple of times to try to verify. I believe that I did but it's definitely a possibility that I have a typo. – R'mani Haulcy Oct 16 '19 at 17:05
• Check first $\nabla.(M_1\nabla H)=(\nabla M_1)(\nabla H)+M_1\nabla^2 H$ (You have 4 typos there eq 11-14). And Abs(D[n[x,y,t],t]) should be Abs[D[n[x,y,t],t]]  (You have 1 typo there eq 5)) – Alex Trounev Oct 16 '19 at 17:24
• I fixed eq 5 by replacing the brackets and I fixed equations 11-14 by replacing the last term by Div[m1[x,y,t] * Grad[h[x,y,t],{x,y,t}],{x,y,t}] or the equivalent. I also used the Laplacian function for the 6 equations that have it instead of computing that by hand. Now, I'm getting the following error: NDSolve: Boundary values may only be specified for one independent variable. Initial values may only be specified at one value of the other independent variable. – R'mani Haulcy Oct 16 '19 at 17:52
• I also updated the code above. – R'mani Haulcy Oct 16 '19 at 17:55

After all the corrections, I also added artificial viscosity (but have not used it yet). The authors of the article did not indicate a solution method. They use average parameters. I used the parameters at point {.5,.5}. Therefore, there is no complete coincidence with Fig. 1 from the article.

(*Parameters*)DAO = 4.32*10^-2; aVis = 0 (*artificial viscosity *);
DH = 8.11*10^-2;
DTAlpha = 6.55*10^-2;
DTBeta = 6.55*10^-2;
Di10 = 6.04*10^-2;
DP = 1.2*10^-1;
lambdaBetaI = 9.51*10^-6;
lambdaN = 8*10^-9;
lambdaA = 8*10^-10;
lambdaTau0 = 8.1*10^-11;
lambdaTau = 1.35*10^-11;
lambdaF = 1.662*10^-3;
lambdaATAlpha = 1.54;
lambdaAABeta0 = 1.793;
lambdaAO = 5*10^-2;
lambdaH = 3*10^-5;
lambdaMF = 2*10^-2;
lambdaMA = 2.3*10^-3;
lambdaM1TBeta = 6*10^-3;
lambdaM1HatTBeta = 6*10^-4;
lambdaTBetaM2 = 1.5*10^-2;
lambdaTBetaM2Hat = 1.5*10^-2;
lambdaTAlphaM1 = 3*10^-2;
lambdaTAlphaM1Hat = 3*10^-2;
lambdai10M2 = 6.67*10^-3;
lambdai10M2Hat = 6.67*10^-3;
lambdaPA = 6.6*10^-8;
lambdaPM2 = 1.32*10^-7;
myTheta = 0.9;
alpha1 = 5;
myBeta = 10;
myGamma = 1;
dABetaI = 9.51;
dABeta0M = 9.51;
dABeta0M = 2*10^-3;
dABeta0MHat = 10^-2;
dTau = 0.277;
dFI = 2.77*10^-3;
dFO = 2.77*10^-4;
dN = 1.9*10^-4;
dNF = 3.4*10^-4;
dNT = 1.7*10^-4;
dNdM = 0.06;
dNdMhat = 0.02;
dA = 1.2*10^-3;
dM1 = 0.015;
dM2 = 0.015;
dM1Hat = 0.015;
dM2Hat = 0.015;
dAO = 0.951;
dH = 58.71;
dTAlpha = 55.45;
dTBeta = 3.33*10^2;
di10 = 16.64;
dP = 1.73;
r0 = 6;
m0 = 5*10^-2;
n0 = 0.14;
mG0 = 0.047;
a0 = 0.14;
kBarABeta0 = 7*10^-3;
kBarNd = 10^-3;
ki10 = 2.5*10^-6;
kTBeta = 2.5*10^-7;
kM = 0.047;
kMHat = 0.047;
kM1 = 0.03;
kM2 = 0.017;
kM1Hat = 0.04;
kM2Hat = 0.007;
kfI = 3.36*10^-10;
kFO = 2.58*10^-11;
kAO = 10^-7;
kP = 6*10^-9;
ktAlpha = 4*10^-5;
fff[t_] := Piecewise[{{r*t/100, t <= 100}, {r, True}}] /. r -> r0;
(*System of Nonlinear PDEs*)
pde = {D[aBetaI[x, y, t],
t] - ((lambdaBetaI*(1 + fff[t]) -
dABetaI*aBetaI[x, y, t])*(n[x, y, t]/n0)) -
aVis Laplacian[aBetaI[x, y, t], {x, y}],
D[aBeta0[x, y, t],
t] - (aBetaI[x, y, t]*
Abs[-dNF*(fI[x, y, t]/(fI[x, y, t] + kfI))*n[x, y, t] -
dNT*(tAlpha[x, y,

t]/(tAlpha[x, y, t] + ktAlpha))*(1/(1 +
myGamma*(i10[x, y, t]/ki10)))*n[x, y, t]] +
lambdaN*(n[x, y, t]/n0) +
lambdaA*(a[x, y, t]/
a0) - (dABeta0MHat*(m1Hat[x, y, t] +
myTheta*m2Hat[x, y, t]) +
dABeta0M*(m1[x, y, t] + myTheta*m2[x, y, t]))*(aBeta0[x, y,
t]/(aBeta0[x, y, t] + kBarABeta0))) -
aVis Laplacian[aBeta0[x, y, t], {x, y}],
D[myTau[x, y, t],
t] - ((lambdaTau0 + lambdaTau*fff[t] -
dTau*myTau[x, y, t])*(n[x, y, t]/n0)),
D[fI[x, y, t],
t] - (lambdaF*myTau[x, y, t] - dFI*fI[x, y, t])*(n[x, y, t]/n0) -
aVis Laplacian[fI[x, y, t], {x, y}],
D[fO[x, y, t],
t] - (fI[x, y, t]*
Abs[-dNF*(fI[x, y, t]/(fI[x, y, t] + kfI))*n[x, y, t] -
dNT*(tAlpha[x, y,
t]/(tAlpha[x, y, t] + ktAlpha))*(1/(1 +
myGamma*(i10[x, y, t]/ki10)))*n[x, y, t]] -
dFO*fO[x, y, t]) - aVis Laplacian[fO[x, y, t], {x, y}],
D[n[x, y, t],
t] - (-dNF*(fI[x, y, t]/(fI[x, y, t] + kfI))*n[x, y, t] -
dNT*(tAlpha[x, y,
t]/(tAlpha[x, y, t] + ktAlpha))*(1/(1 +
myGamma*(i10[x, y, t]/ki10)))*n[x, y, t]),
D[a[x, y, t],
t] - (lambdaAABeta0*aBeta0[x, y, t] +
lambdaATAlpha*tAlpha[x, y, t] - dA*a[x, y, t]) -
aVis Laplacian[a[x, y, t], {x, y}],
D[nD[x, y, t],
t] - (dNF*(fI[x, y, t]/(fI[x, y, t] + kfI))*n[x, y, t] +
dNT*(tAlpha[x, y,
t]/(tAlpha[x, y, t] + ktAlpha))*(1/(1 +
myGamma*(i10[x, y, t]/ki10)))*n[x, y, t] -
dNdM*(m1[x, y, t] +
m2[x, y, t])*(nD[x, y, t]/(nD[x, y, t] + kBarNd)) -
dNdMhat*(m1Hat[x, y, t] +
m2Hat[x, y, t])*(nD[x, y, t]/(nD[x, y, t] + kBarNd))) -
aVis Laplacian[nD[x, y, t], {x, y}],
D[aO[x, y, t],
t] - ((lambdaAO*aBeta0[x, y, t] - dAO*aO[x, y, t] +
DAO*(D[aO[x, y, t], {x, 2}] + D[aO[x, y, t], {y, 2}]))),
D[h[x, y, t],
t] - (lambdaH*nD[x, y, t] - dH*h[x, y, t] +
DH*(D[h[x, y, t], {x, 2}] + D[h[x, y, t], {y, 2}])),
D[m1[x, y, t],
t] - (mG0*(lambdaMF*(fO[x, y, t]/(fO[x, y, t] + kFO)) +
lambdaMA*(aO[x, y,
t]/(aO[x, y, t] + kAO)))*((myBeta*(tAlpha[x, y,

t]/(tAlpha[x, y, t] + ktAlpha)))/(myBeta*(tAlpha[x, y,
t]/(tAlpha[x, y, t] + ktAlpha)) + (i10[x, y,
t]/(i10[x, y, t] + ki10)))) -
lambdaM1TBeta*(tBeta[x, y, t]/(tBeta[x, y, t] + kTBeta))*
m1[x, y, t] -
dM1*(D[m1[x, y, t], x] D[h[x, y, t], x] +
D[m1[x, y, t], y] D[h[x, y, t], y] +
m1[x, y,
t]*(D[h[x, y, t], {x, 2}] + D[h[x, y, t], {y, 2}]))),
D[m2[x, y, t],
t] - (mG0*(lambdaMF*(fO[x, y, t]/(fO[x, y, t] + kFO)) +
lambdaMA*(aO[x, y, t]/(aO[x, y, t] + kAO)))*((i10[x, y,
t]/(i10[x, y, t] +
ki10))/(myBeta*(tAlpha[x, y,
t]/(tAlpha[x, y, t] + ktAlpha)) + (i10[x, y,
t]/(i10[x, y, t] + ki10)))) +
lambdaM1TBeta*(tBeta[x, y, t]/(tBeta[x, y, t] + kTBeta))*
m1[x, y, t] -
dM2*(D[m2[x, y, t], x] D[h[x, y, t], x] +
D[m2[x, y, t], y] D[h[x, y, t], y] +
m2[x, y,
t]*(D[h[x, y, t], {x, 2}] + D[h[x, y, t], {y, 2}]))),
D[m1Hat[x, y, t],
t] - ((alpha1*(p[x, y, t]/(p[x, y, t] + kP)))*(m0 - (m1Hat[x, y,
t] + m2Hat[x, y,

t]))*((myBeta*(tAlpha[x, y,
t]/(tAlpha[x, y, t] + ktAlpha)))/(myBeta*(tAlpha[x, y,
t]/(tAlpha[x, y, t] + ktAlpha)) + (i10[x, y,
t]/(i10[x, y, t] + ki10)))) -
lambdaM1HatTBeta*(tBeta[x, y, t]/(tBeta[x, y, t] + kTBeta))*
m1Hat[x, y, t] -
dM1Hat*(D[m1Hat[x, y, t], x] D[aO[x, y, t], x] +
D[m1Hat[x, y, t], y] D[aO[x, y, t], y] +
m1Hat[x, y,
t]*(D[aO[x, y, t], {x, 2}] + D[aO[x, y, t], {y, 2}]))),
D[m2Hat[x, y, t],
t] - ((alpha1*(p[x, y, t]/(p[x, y, t] + kP)))*(m0 - (m1Hat[x, y,
t] + m2Hat[x, y, t]))*((i10[x, y,
t]/(i10[x, y, t] +
ki10))/(myBeta*(tAlpha[x, y,
t]/(tAlpha[x, y, t] + ktAlpha)) + (i10[x, y,
t]/(i10[x, y, t] + ki10)))) +
lambdaM1HatTBeta*(tBeta[x, y, t]/(tBeta[x, y, t] + kTBeta))*
m1Hat[x, y, t] -
dM2Hat*(D[m2Hat[x, y, t], x] D[aO[x, y, t], x] +
D[m2Hat[x, y, t], y] D[aO[x, y, t], y] +
m2Hat[x, y,
t]*(D[aO[x, y, t], {x, 2}] + D[aO[x, y, t], {y, 2}]))),
D[tBeta[x, y, t],
t] - (lambdaTBetaM2*m2[x, y, t] +
lambdaTBetaM2Hat*m2Hat[x, y, t] - dTBeta*tBeta[x, y, t] +
DTBeta*(D[tBeta[x, y, t], {x, 2}] + D[tBeta[x, y, t], {y, 2}])),
D[tAlpha[x, y, t],
t] - (lambdaTAlphaM1*m1[x, y, t] +
lambdaTAlphaM1Hat*m1Hat[x, y, t] - dTAlpha*tAlpha[x, y, t] +
DTAlpha*(D[tAlpha[x, y, t], {x, 2}] +
D[tAlpha[x, y, t], {y, 2}])),
D[i10[x, y, t],
t] - (lambdai10M2*m2[x, y, t] + lambdai10M2Hat*m2Hat[x, y, t] -
di10*i10[x, y, t] +
Di10*(D[i10[x, y, t], {x, 2}] + D[i10[x, y, t], {y, 2}])),
D[p[x, y, t],
t] - (lambdaPA*a[x, y, t] + lambdaPM2*m2[x, y, t] -
dP*p[x, y, t] +
DP*(D[p[x, y, t], {x, 2}] + D[p[x, y, t], {y, 2}]))};

(*Periodic Boundary Conditions*)
bc = {aO[0, y, t] == aO[1, y, t], aO[x, 0, t] == aO[x, 1, t],
h[0, y, t] == h[1, y, t], h[x, 0, t] == h[x, 1, t],
tBeta[0, y, t] == tBeta[1, y, t], tBeta[x, 0, t] == tBeta[x, 1, t],
i10[0, y, t] == i10[1, y, t], i10[x, 0, t] == i10[x, 1, t],
tAlpha[0, y, t] == tAlpha[1, y, t],
tAlpha[x, 0, t] == tAlpha[x, 1, t], p[0, y, t] == p[1, y, t],
p[x, 0, t] == p[x, 1, t]};

(*Initial Conditions*)
ic = {aBetaI[x, y, 0] == 10^(-6), aBeta0[x, y, 0] == 10^(-8),
myTau[x, y, 0] == 1.37*10^(-10), fI[x, y, 0] == 3.36*10^(-10),
fO[x, y, 0] == 3.36*10^(-11), n[x, y, 0] == 0.14,
a[x, y, 0] == 0.14, m1[x, y, 0] == 0.02, m2[x, y, 0] == 0.02,
m1Hat[x, y, 0] == 0, m2Hat[x, y, 0] == 0, nD[x, y, 0] == 0,
h[x, y, 0] == 1.3*10^(-11), tBeta[x, y, 0] == 10^(-6),
tAlpha[x, y, 0] == 2*10^(-5), i10[x, y, 0] == 10^(-5),
p[x, y, 0] == 5*10^(-9), aO[x, y, 0] == 0};

eqns = Flatten@{Table[pde[[i]] == 0, {i, 1, Length[pde]}], bc,
ic}; var = {aBetaI, aBeta0, myTau, fI, fO, n, a, m1, m2, m1Hat,
m2Hat, nD, h, tBeta, tAlpha, i10, p, aO};

Monitor[sols =
NDSolve[eqns, var, {x, 0, 1}, {y, 0, 1}, {t, 0, 3650},
Method -> {"IndexReduction" -> Automatic,
"EquationSimplification" -> "Residual",
"PDEDiscretization" -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid",
"MinPoints" -> 41, "MaxPoints" -> 41,
"DifferenceOrder" -> 2}}},
EvaluationMonitor :> (monitor =
Row[{"t=", CForm[t], " csol=", CForm[n[.5, .5, t]]}])],
monitor];


Solution of the system of equations at point {x,y}={.5,.5}

point[i_] :=
With[{max =
Max[Table[var[[i]][.5, .5, t] /. sols[[1, i]], {t, 0, 3650, 50}]],
min = Min[
Table[var[[i]][.5, .5, t] /. sols[[1, i]], {t, 0, 3650, 50}]]},
Range[min, max, (-min + max)/4]]

Table[Plot[var[[i]][.5, .5, t] /. sols[[1, i]], {t, 0, 3650},
PlotLabel -> var[[i]], PlotRange -> All,
Ticks -> {{0, 2000, 3650}, point[i]}], {i, 1, Length[var]}]


Solution of the system of equations at point t=3650

Table[With[{max =
Max[Table[var[[i]][.5, .5, t] /. sols[[1, i]], {t, 0, 3650, 50}]],
min = Min[
Table[var[[i]][.5, .5, t] /. sols[[1, i]], {t, 0, 3650, 50}]]},
Plot3D[var[[i]][x, y, 3650] /. sols[[1, i]], {x, 0, 1}, {y, 0, 1},
PlotLabel -> var[[i]], PlotRange -> {min, max}, Mesh -> None,
ColorFunction -> Hue]], {i, 1, Length[var]}]


• Thank you so much for all of your help. I need the plots to match, so I will try to use the code you've provided to get matching plots. – R'mani Haulcy Oct 18 '19 at 13:26
• @R'maniHaulcy You're welcome! – Alex Trounev Oct 18 '19 at 14:34