Параболическое + эллиптическое уравнение в Wolfram Mathematica

Недавно я писал про решение системы параболического и эллиптического уравнений в MATLAB (см. здесь). Выяснилось, что MATLAB не справился с задачей решения такой системы, когда для функций ставятся краевые условия 3-го рода.
Попробуем решить такую же систему в Wolfram Mathematica.
Если записать такую систему непосредственно в функции NDSolve, то Mathematica её не решит. Однако для решения систем такого типа (когда система дифференциальных уравнений в частных производных методом прямых сводится к обыкновенным дифференциальным уравнениям вместе с алгебраическими уравнениями) можно использовать подход, описанный в статье Numerical Solution of Differential-Algebraic Equations документации Wolfram Mathematica (см. раздел Combined Elliptic-Parabolic PDE in 1D).

Исходный код (это немного изменённый пример, который приведён в той статье):

ll = 1;

Clear[u, v, nx, eqn1, eqn2, ic1, ic2, vars,
  eqn1, eqn2, ic1, ic2, vars, sol, PDEsol, order, p, npts];
npts = 100; order = 12; p = 10;
nx = Range[0, ll, 1/(npts - 1)];
{ddx, d2dx2} =
  Map[NDSolve`FiniteDifferenceDerivative[Derivative[#], nx,
      "DifferenceOrder" -> order]["DifferentiationMatrix"] &, {1, 2}];
u = Array[u$[#][t] &, npts]; v = Array[v$[#][t] &, npts];
eqn1 = Thread[D[u, t] == d2dx2.u + 2];
eqn2 = Thread[0 == d2dx2.v + u];

eqn1[[1]] = -ddx[[1]].u + u[[1]] == 2*t;
eqn1[[-1]] = ddx[[-1]].u + u[[-1]] == 2*t;

eqn2[[1]] = -ddx[[1]].v + v[[1]] == 0;
eqn2[[-1]] = ddx[[-1]].v + v[[-1]] == 0;

ic1 = Thread[u == 0] /. t -> 0;
ic2 = Thread[v == 0] /. t -> 0;

vars = Flatten[{u, v}] /. x__[t] :> x;
PDEsol = First[NDSolve[{eqn1, eqn2, ic1, ic2}, vars, {t, 0, 3}]];

time = Flatten@First[vars[[1]] /. PDEsol];
{usol, vsol} =
  Flatten[
     Table[Transpose[{nx, ConstantArray[ti, npts], # /. PDEsol}] /.
       t -> ti,
      {ti, time[[1]], time[[2]], 0.1}], 1] & /@ {u, v};
GraphicsRow[ListPlot3D[#, PlotRange -> All] & /@ {usol, vsol},
 ImageSize -> 500]

Wolfram Mathematica показала отличные результаты. Вот графики приближённых решений.

График приближённого решения u(x,t)

График приближённого решения v(x,t)

Отметим, что для того, чтобы начальное условие могло зависеть от x и чтобы в уравнении могла стоять функция от x, этот пример нужно ещё как-то подправить.
Вот исходный код для решения той же системы, но со сдвинутым на 1 временем (в этом случае начальное условие для v зависит от x).

ll = 1;

Clear[u, v, nx, eqn1, eqn2, ic1, ic2, vars,
  eqn1, eqn2, ic1, ic2, vars, sol, PDEsol, order, p, npts];
npts = 100; order = 12; p = 10;
nx = Range[0, ll, 1/(npts - 1)];
{ddx, d2dx2} =
  Map[NDSolve`FiniteDifferenceDerivative[Derivative[#], nx,
      "DifferenceOrder" -> order]["DifferentiationMatrix"] &, {1, 2}];
u = Array[u$[#][t] &, npts]; v = Array[v$[#][t] &, npts];
eqn1 = Thread[D[u, t] == d2dx2.u + 2];
eqn2 = Thread[0 == d2dx2.v + u];

eqn1[[1]] = -ddx[[1]].u + u[[1]] == 2*(t + 1);
eqn1[[-1]] = ddx[[-1]].u + u[[-1]] == 2*(t + 1);

eqn2[[1]] = -ddx[[1]].v + v[[1]] == 0;
eqn2[[-1]] = ddx[[-1]].v + v[[-1]] == 0;

ic1 = Thread[(u /. t -> 0) == 2];
ic2 = Thread[(v /. t -> 0) == (nx /. x_ :> (1 + x - x^2))];

vars = Flatten[{u, v}] /. x__[t] :> x;
PDEsol = First[NDSolve[{eqn1, eqn2, ic1, ic2}, vars, {t, 0, 3}]];

time = Flatten@First[vars[[1]] /. PDEsol];
{usol, vsol} =
  Flatten[
     Table[Transpose[{nx, ConstantArray[ti, npts], # /. PDEsol}] /.
       t -> ti,
      {ti, time[[1]], time[[2]], 0.1}], 1] & /@ {u, v};
GraphicsRow[ListPlot3D[#, PlotRange -> All] & /@ {usol, vsol},
 ImageSize -> 500]

uu[tt_] := (u /. PDEsol) /. t -> tt;

ures[xx_, tt_] := ListInterpolation[uu[tt], {0, ll}][xx];

vv[tt_] := (v /. PDEsol) /. t -> tt;

vres[xx_, tt_] := ListInterpolation[vv[tt], {0, ll}][xx];

Plot3D[ures[x, t], {x, 0, ll}, {t, 0, 3}]
Plot3D[vres[x, t], {x, 0, ll}, {t, 0, 3}]