Here's a hacky way to do it, based on inspecting this:
NDSolve`ProcessEquations[{eq1, iv5, bcs}, {u[x, t]}, {x, 0, Xmax}, {t, 0, Tmax}]
There's a couple of places where you see MapThread[rhsFN, data, 1]
, that maps the right-hand side of the first-orderized differential equation onto the state data. Since in this case, the RHS is vectorized, we can override MapThread
and apply the RHS directly with a integration slipped in for a dummy function int[]
. Maybe not the safest way to do this, but I thought it was cool enough to share.
Xmax = 5;
Tmax = 5;
eq1 = D[u[x, t], t] == D[u[x, t], x, x] + (x - int[u[x, t], x, t]) u[x, t]; (* N.B. *)
iv5 = {u[x, 0] == 2/(Sqrt[Pi]*Exp[x^2])};
bcs = {u[0, t] == 2/Sqrt[Pi], u[Xmax, t] == 0};
Block[{int, xx},
int[u_, x_, t_ /; t == 0] = (* IC - fools ProcessEquations, thinks int[] a good num.fn. *)
NIntegrate[2/(Sqrt[Pi]*Exp[x^2]), {x, 0, Xmax}];
int[u_?VectorQ, x_?VectorQ, t_?NumericQ] :=
Integrate[Interpolation[Transpose@{x, x*u}][xx], xx] /. xx -> Xmax;
Internal`InheritedBlock[{MapThread},
Unprotect[MapThread];
MapThread[f_, data_, 1] /; ! FreeQ[f, int] := f @@ data;
Protect[MapThread];
s10 = NDSolve[{eq1, iv5, bcs}, {u[x, t]}, {x, 0, Xmax}, {t, 0, Tmax}];
]];
Plot3D[u[x, t] /. s10, {x, 0, 5}, {t, 0, 5}]