# How to plot a system of coupled PDEs using NDSolve [closed]

I am trying to solve the the following system of couples PDEs using NDSolve. However I am struggling to plot the solutions. I am trying to plot |u1[t,x]|^2 (i.e. the square of the absolute value) but whenever I try I get a blank 3D box as opposed to any actual results. How could I get the results to show in 3D? I have attached the code below.

I have tried the following line:

Plot3D[u1[t, x] /. sol, {t, 0, 5}, {x, -a, a}] but it generates an blank box.


eq = {

I*D[u1[t, x], {t, 1}] + D[u1[t, x], {x, 2}] + u6[t,x]+u2[t,x]  == 0,
I*D[u2[t, x], {t, 1}] + D[u2[t, x], {x, 2}] + u1[t,x]+u3[t,x]  == 0,
I*D[u3[t, x], {t, 1}] + D[u3[t, x], {x, 2}] + u2[t,x]+u4[t,x]  == 0,
I*D[u4[t, x], {t, 1}] + D[u4[t, x], {x, 2}] + u3[t,x]+u5[t,x]  == 0,
I*D[u5[t, x], {t, 1}] + D[u5[t, x], {x, 2}] + u4[t,x]+u6[t,x]  == 0,
I*D[u6[t, x], {t, 1}] + D[u6[t, x], {x, 2}] + u5[t,x]+u1[t,x]  == 0,

u1[0, x] == f1,
u2[0, x] == f2,
u3[0, x] == f3,
u4[0, x] == f4,
u5[0, x] == f5,
u6[0, x] == f6,

u1[t, -a] == u1[t, a],
u2[t, -a] == u2[t, a],
u3[t, -a] == u3[t, a],
u4[t, -a] == u4[t, a],
u5[t, -a] == u5[t, a],
u6[t, -a] == u6[t, a]

};

sol = NDSolve[eq,{u1[t,x],u2[t,x],u3[t,x],u4[t,x],u5[t,x],u6[t,x]}, {t, 0, 5}, {x, -a, a}, MaxStepSize -> 0.07];


• This is not the complete code. What values did you use for a, f1, f2, f3, f4, f5, f6? – Bob Hanlon Jul 24 at 12:29
• a has value 10 and f1 = 0, f2 = 0,f3 = 0, f4 = 0, f5 = 0, f6 = 0 for example, I am just trying to display something as opposed to getting a blank box. – Dman Jul 24 at 12:32
• Edit your question to include the values given in your comment. Plot3D has the attribute HoldAll. Use Evaluate, i.e., Plot3D[Evaluate[u1[t, x] /. sol], {t, 0, 5}, {x, -a, a}] – Bob Hanlon Jul 24 at 12:37

We can't get something from zero initial data and periodic boundary condition. If we put some initial data, for instance for u1, then we have (see attentively how we use sol and Abs for complex solution):

eq = {I*D[u1[t, x], {t, 1}] + D[u1[t, x], {x, 2}] + u6[t, x] +
u2[t, x] == 0,
I*D[u2[t, x], {t, 1}] + D[u2[t, x], {x, 2}] + u1[t, x] +
u3[t, x] == 0,
I*D[u3[t, x], {t, 1}] + D[u3[t, x], {x, 2}] + u2[t, x] +
u4[t, x] == 0,
I*D[u4[t, x], {t, 1}] + D[u4[t, x], {x, 2}] + u3[t, x] +
u5[t, x] == 0,
I*D[u5[t, x], {t, 1}] + D[u5[t, x], {x, 2}] + u4[t, x] +
u6[t, x] == 0,
I*D[u6[t, x], {t, 1}] + D[u6[t, x], {x, 2}] + u5[t, x] +
u1[t, x] == 0, u1[0, x] == Sin[2 Pi x/a], u2[0, x] == f2,
u3[0, x] == f3, u4[0, x] == f4, u5[0, x] == f5, u6[0, x] == f6,
u1[t, -a] == u1[t, a], u2[t, -a] == u2[t, a],
u3[t, -a] == u3[t, a], u4[t, -a] == u4[t, a],
u5[t, -a] == u5[t, a], u6[t, -a] == u6[t, a]};
par = {a = 10, f1 = 1, f2 = 0, f3 = 0, f4 = 0, f5 = 0,
f6 = 0}; var = {u1, u2, u3, u4, u5, u6};
sol = NDSolve[eq, var, {t, 0, 5}, {x, -a, a},
Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid",
"MinPoints" -> 100, "MaxPoints" -> 100,
"DifferenceOrder" -> "Pseudospectral"}}]

Table[Plot3D[
Evaluate[Abs[var[[i]][t, x]] /. First[sol]], {t, 0, 5}, {x, -a, a},
Mesh -> None, ColorFunction -> "Rainbow", AxesLabel -> Automatic,
PlotLabel -> var[[i]]], {i, Length[var]}]