| Вот целиком код:
 c1 = 0.1; c2 = 90; c3 = 10; c4 = 0.1; eps = 0.025; k1 = 5; k2 = 20; T= 10;
 
 n[g_] = c1 + c2 g/(c3 + g);
 v[g_] = c4/(c4 + g);
 l[a_, b_] = (a - b)/(a + b);
 
 eqn1 = D[Pr[x, t], t] == D[Pr[x, t], x] + n[G[x, t]];
 eqn2 = D[Pl[x, t], t] == -D[Pl[x, t], x] + n[G[x, t]];
 eqn3 = D[Mr[x, t], t] == D[v[G[x, t]] Mr[x, t], x] + n[G[x, t]];
 eqn4 = D[Ml[x, t], t] == -D[v[G[x, t]] Ml[x, t], x] + n[G[x, t]];
 
 eqn6 = eps D[Gr[x, t], t] == -k1 Gr[x, t] + k2 Nr[x, t] Gs[x, t] - D[Gr[x, t], x];
 eqn7 = eps D[Gl[x, t], t] == -k1 Gl[x, t] + k2 Nl[x, t] Gs[x, t] + D[Gl[x, t], x];
 eqn8 = eps D[Gst[x, t], t] == k1 (Gl[x, t] + Gr[x, t]) - k2 (Nl[x, t] + Nr[x, t]) Gst[x, t];
 
 G[x_, t_] = 0.5 (1 + x^2)^(-k1);
 
 
 sol1 = NDSolve[{eqn1, Pr[1, t] == 0, Pr[x, 0] == 1 - x}, Pr[x, t], {x, -1, 1}, {t, 0, T}, PrecisionGoal -> 3]
 
 sol2 = NDSolve[{eqn2, Pl[-1, t] == 0, Pl[x, 0] == x + 1}, Pl[x, t], {x, -1, 1}, {t, 0, T}, PrecisionGoal -> 3]
 
 sol3 = NDSolve[{eqn3, Mr[1, t] == 0, Mr[x, 0] == 2 - 2 x},Mr[x, t], {x, -1, 1}, {t, 0, T}, PrecisionGoal -> 2]
 
 sol4 = NDSolve[{eqn4, Ml[-1, t] == 0, Ml[x, 0] == 2 x + 2}, Ml[x, t], {x, -1, 1}, {t, 0, T}, PrecisionGoal -> 2]
 
 rr[x_, t_] = Evaluate[Mr[x, t] /. sol3];
 rrr[x_, t_] = Evaluate[Pr[x, t] /. sol1];
 ll[x_, t_] = Evaluate[Ml[x, t] /. sol4];
 lll[x_, t_] = Evaluate[Pl[x, t] /. sol2];
 Nl[x_, t_] = Integrate[(ll[y, t][[1]] - lll[y, t][[1]]), {y, -1, x}];
 Nr[x_, t_] = Integrate[(rr[y, t][[1]] - rrr[y, t][[1]]), {y, x, 1}];
 
 sol678 = NDSolve[{eqn6, eqn7, eqn8,
 Gr[-1, t] == 0, Gl[1, t] == 0, Gst[-1, t] == 1/64,
 Gr[x, 0] == 0, Gl[x, 0] == 0, Gst[x, 0] == 0.5 (x^2 + 1)^(-k1)},
 {Gr[x, t], Gl[x, t], Gst[x, t]}, {x, -1, 1}, {t, 0, T}]
 
 |