sol61115 = NDSolve[{D[u1[x, t], t] == 1*D[u1[x, t], {x, 2}] - 1*u1[x, t]^2/(u1[x, t]^2 + 1^2) - 1*u1[x, t]^2/(u1[x, t]^2 + 1^2) + 1*1*Exp[((1*1^2/(1^2 + 1^2) + 1 - 1))/1] + 1*1*1*u1[x, t]^2*u2[x, t]^2* u3[x, t]/((u1[x, t]^2 + 1^2)*(u2[x, t]^2 + 1^2)) + 1*1*1*u1[x, t]^2*u2[x, t]^2* u4[x, t]/((u1[x, t]^2 + 1^2)*(u2[x, t]^2 + 1^2)),(*1式*) D[u2[x, t], t] == 1*D[u2[x, t], {x, 2}] - (1*1/(1 + 1))*u1[x, t]/(u1[x, t] + 1) - 1*u2[x, t],(*2式*)(*边界条件*) 1*D[u3[x, t], t] == 1^2/(1^2 + u1[x, t]^2) - u3[x, t], 1*D[u4[x, t], t] == 1^2/(1^2 + u1[x, t]^2) - u4[x, t], D[u1[0, t], t] == 0, D[u1[1, t], t] == 0, D[u2[0, t], t] == 0, D[u1[1, t], t] == 0},(*初值条件*){u1, u2, u3, u4}, {x, 0, 1}, {t, 0, 1}]

Plot3D[Evaluate[{u1[x, t], u2[x, t], u3[x, t], u4[x, t]} /. sol61115], {t, 0, 1}, {x, 0, 1}, PlotRange -> All, AxesLabel -> {x, t, u}, PlotStyle -> {{Opacity[.3], Blue}, {Opacity[.3], Red}, {Opacity[.3], Black}, {Opacity[.3], Pink}}, MeshFunctions -> {#3 &}]
这个是绘图
自己是刚从MATLAB转过来的小白，不太懂，跪谢各位大佬

稍微修改了下错误减少了，但是提示出现了NDSolve::femcnsd

sol61115 = NDSolve[{D[u1[x, t], t] == 1*D[u1[x, t], {x, 2}] - 1*u1[x, t]^2/(u1[x, t]^2 + 1^2) - 1*u1[x, t]^2/(u1[x, t]^2 + 1^2) + 1*1*(1 + exp[((1*1^2/(1^2 + 1^2) + 1 - 1))/1]) + 1*1*1*u1[x, t]^2*u2[x, t]^2* u3[x, t]/((u1[x, t]^2 + 1^2)*(u2[x, t]^2 + 1^2)) + 1*1*1*u1[x, t]^2*u2[x, t]^2* u4[x, t]/((u1[x, t]^2 + 1^2)*(u2[x, t]^2 + 1^2)),(*1式*) D[u2[x, t], t] == 1*D[u2[x, t], {x, 2}] - (1*1/(1 + 1))*u1[x, t]/(u1[x, t] + 1) - 1*u2[x, t],(*2式*)(*边界条件*) 1*D[u3[x, t], t] == 1^2/(1^2 + u1[x, t]^2) - u3[x, t](*3式*), 1*D[u4[x, t], t] == 1^2/(1^2 + u1[x, t]^2) - u4[x, t](*4式*), D[u1[0, t], t] == 0, D[u1[1, t], t] == 0, D[u2[0, t], t] == 0, D[u2[1, t], t] == 0},(*初值条件*){u1, u2, u3, u4}, {x, 0, 1}, {t, 0, 1}]
再改了一下，NDsolve好像解不开..

自己加了t=0的边界条件，但是还是解不开，一直说true不能做变量？我代码好像也没true在里面

sol61115 = NDSolve[{D[u1[x, t], t] == 1*D[u1[x, t], {x, 2}] + 0.1 - 1*u1[x, t]^2/(u1[x, t]^2 + 1^2) - 1*u1[x, t]^2/(u1[x, t]^2 + 1^2) + 1*1*(1 + exp[((1*1^2/(1^2 + 1^2) + 1 - 1))/1]) + 1*1*1*u1[x, t]^2*u2[x, t]^2* u3[x, t]/((u1[x, t]^2 + 1^2)*(u2[x, t]^2 + 1^2)) + 1*1*1*u1[x, t]^2*u2[x, t]^2* u4[x, t]/((u1[x, t]^2 + 1^2)*(u2[x, t]^2 + 1^2)),(*1式*) D[u2[x, t], t] == 1*D[u2[x, t], {x, 2}] - (1*1/(1 + 1))*u1[x, t]/(u1[x, t] + 1) - 1*u2[x, t],(*2式*) 1*D[u3[x, t], t] == 1^2/(1^2 + u1[x, t]^2) - u3[x, t](*3式*), 1*D[u4[x, t], t] == 1^2/(1^2 + u1[x, t]^2) - u4[x, t](*4式*), u1[x, 0] == 0, u2[x, 0] == 0, u3[x, 0] == 0, u4[x, 0] == 0, D[u1[0, t], t] == 0, D[u1[1, t], t] == 0, D[u2[0, t], t] == 0, D[u2[1, t], t] == 0},(*边界条件*){u1, u2, u3, u4}, {x, 0, 1}, {t, 0, 1}]Plot3D[Evaluate[{u1[x, t], u2[x, t], u3[x, t], u4[x, t]} /. sol61115], {t, 0, 1}, {x, 0, 1}, PlotRange -> All, AxesLabel -> {x, t, u}, PlotStyle -> {{Opacity[.3], Blue}, {Opacity[.3], Red}, {Opacity[.3], Black}, {Opacity[.3], Pink}}, MeshFunctions -> {#3 &}]

sol61115 = NDSolve[{D[u1[x, t], t] == 1*D[u1[x, t], {x, 2}] + 0.1 - 1*u1[x, t]^2/(u1[x, t]^2 + 1^2) - 1*u1[x, t]^2/(u1[x, t]^2 + 1^2) + 1*1*(1 + Exp[((1*1^2/(1^2 + 1^2) + 1 - 1))/1]) + 1*1*1*u1[x, t]^2*u2[x, t]^2* u3[x, t]/((u1[x, t]^2 + 1^2)*(u2[x, t]^2 + 1^2)) + 1*1*1*u1[x, t]^2*u2[x, t]^2* u4[x, t]/((u1[x, t]^2 + 1^2)*(u2[x, t]^2 + 1^2)),(*1式*) D[u2[x, t], t] == 1*D[u2[x, t], {x, 2}] - (1*1/(1 + 1))*u1[x, t]/(u1[x, t] + 1) - 1*u2[x, t],(*2式*) 1*D[u3[x, t], t] == 1^2/(1^2 + u1[x, t]^2) - u3[x, t](*3式*), 1*D[u4[x, t], t] == 1^2/(1^2 + u1[x, t]^2) - u4[x, t](*4式*), u1[x, 0] == 0, u2[x, 0] == 0, u3[x, 0] == 0, u4[x, 0] == 0, Derivative[1, 0][u1][0, t] == 0, Derivative[1, 0][u1][1, t] == 0, Derivative[1, 0][u2][0, t] == 0, Derivative[1, 0][u2][1, t] == 0},(*边界条件*){u1, u2, u3, u4}, {x, 0, 1}, {t, 0, 1}]Clear@DerivativePlot3D[Evaluate[{u1[x, t], u2[x, t], u3[x, t], u4[x, t]} /. sol61115], {t, 0, 1}, {x, 0, 1}, PlotRange -> All, AxesLabel -> {x, t, u}, PlotStyle -> {{Opacity[.3], Blue}, {Opacity[.3], Red}, {Opacity[.3], Black}, {Opacity[.3], Pink}}, MeshFunctions -> {#3 &}]
绘制结果已经出了，绘制的很成功，非常感谢吧主大佬

真的神器