We can solve this problem in two steps. First we compute fractional Laplacian $v=(-\Delta)^\frac{s}{2}u$, and then we solve fractional thin film equation
$\partial_tu=\partial_x(u^k\partial_x v)$. Code for the step one is based on the colocation method with using Euler wavelets (see my previous answer here and here). Let define exact solution to test code
c[n_, a_] := a 2^(a - 1) Gamma[(a + n)/2]/(Pi^(n/2) Gamma[1 - a/2]);
lap[n_, a_, x_, u_] :=
c[n, a] Integrate[(u[x] - u[y])/Abs[x - y]^(n + a), y];
ue[a_, x_] := (1 + x^2)^(-(1 - a)/2);
uc[x_, t_] := ue[s, x] (c0 kp (1/(c0 kp) - t))^(-1/kp);
vc[x_, t_] := (c0 kp (1/(c0 kp) - t))^(-1/kp)
2^s Gamma[(1 + s)/2]/Gamma[(1 - s)/2] (1 + x^2)^(-(1 + s)/2);
lape[a_, x_, t_] := (c0 kp (1/(c0 kp) - t))^(-1/kp)
2^a Gamma[(1 + a)/2]/Gamma[(1 - a)/2] (1 + x^2)^(-(1 + a)/2);
f[x_, t_] :=
1/Gamma[(1 - s)/2] (1 - c0 kp t)^(-1 - 1/kp) (1 + x^2)^(-(5/2) - s/
2) (-c0 (1 + x^2)^(2 + s) Gamma[(1 - s)/2] +
2^s kp (-1 + s) (1 + s) (-1 +
c0 kp t) x^2 ((1 - c0 kp t)^(-1/kp) (1 + x^2)^(1/2 (-1 + s)))^
kp Gamma[(1 + s)/2] -
2^s (1 + s) (3 + s) (-1 +
c0 kp t) x^2 ((1 - c0 kp t)^(-1/kp) (1 + x^2)^(1/2 (-1 + s)))^
kp Gamma[(1 + s)/2] +
2^s (1 + s) (-1 + c0 kp t) (1 +
x^2) ((1 - c0 kp t)^(-1/kp) (1 + x^2)^(1/2 (-1 + s)))^
kp Gamma[(1 + s)/2]);
Here kp=k
and c0
some parameters defined below. On these lines we define wavelets and function to be solved
UE[m_, t_] := EulerE[m, t]
psi[k_, n_, m_, t_] :=
Piecewise[{{2^(k/2) Sqrt[2/Pi] UE[m, 2^k t - 2 n + 1], (n - 1)/
2^(k - 1) <= t < n/2^(k - 1)}, {0, True}}]
PsiE[k_, M_, t_] :=
Flatten[Table[psi[k, n, m, t], {n, 1, 2^(k - 1)}, {m, 0, M - 1}]]
k0 = 2; M0 = 3; x0 = 0; x1 = 1; nn =
Total[With[{k = k0, M = M0},
Flatten[Table[1, {n, 1, 2^(k - 1)}, {m, 0, M - 1}]]]]; dx = (x1 -
x0)/(nn); xl = Table[x0 + l*dx, {l, 0, nn}]; xcol =
tcol = Table[(xl[[l - 1]] + xl[[l]])/2, {l, 2, nn + 1}]; Psijk =
With[{k = k0, M = M0}, PsiE[k, M, t1]]; Intx1 =
With[{k = k0, M = M0}, Integrate[PsiE[k, M, t1], t1]];
Intx2 = Integrate[Intx1, t1];
intx1[y_] := Intx1 /. t1 -> (y - x0)/(x1 - x0);
intx2[y_] := Intx2 /. t1 -> (y - x0)/(x1 - x0);
Psi[y_] := Psijk /. t1 -> (y - x0)/(x1 - x0);
varu[t_] := Table[u[i][t], {i, nn}];
varv[t_] := Table[v[i][t], {i, nn}];
v[x_, t_] := varv[t] . intx2[x] + x v0[t] + v01[t];
v1[x_, t_] := varv[t] . intx1[x] + v0[t];
v2[x_, t_] := varv[t] . Psi[x];
u[x_, t_] := varu[t] . intx1[x] + u01[t];
u1[x_, t_] := varu[t] . Psi[x];
Code to compute fractional Laplacian
s = 9/10; kp = 3; c0 = 1/6; int =
Table[Table[
NIntegrate[(intx1[tcol[[i]]][[j]] -
intx1[y][[j]])/(tcol[[i]] - y)^(1 + s), {y, x0, tcol[[i]]},
Method -> "PrincipalValue", Exclusions -> tcol[[i]] - y == 0,
AccuracyGoal -> 6, PrecisionGoal -> 6], {j, nn}], {i, nn}] //
Quiet;
int1 = Table[
Table[NIntegrate[(intx1[tcol[[i]]][[j]] -
intx1[y][[j]])/(y - tcol[[i]])^(1 + s), {y, tcol[[i]], x1},
Method -> "PrincipalValue", Exclusions -> y - tcol[[i]] == 0,
AccuracyGoal -> 6, PrecisionGoal -> 6], {j, nn}], {i, nn}] //
Quiet;
int0 = Table[
Table[NIntegrate[(intx1[tcol[[i]]][[
j]])/(tcol[[i]] - y)^(1 + s), {y, -Infinity, x0},
AccuracyGoal -> 10, PrecisionGoal -> 8], {j, nn}], {i, nn}] //
Quiet;
int2 = Table[
Table[NIntegrate[(intx1[tcol[[i]]][[
j]])/(y - tcol[[i]])^(1 + s), {y, x1, Infinity},
AccuracyGoal -> 10, PrecisionGoal -> 8], {j, nn}], {i, nn}] //
Quiet;
lp = c[1, s] (int + int1 + int0 + int2) // Re;
System of equations, initial and boundary conditions
eq = Join[
Table[D[varu[t], t] . lp[[i]] - D[v[xcol[[i]], t], t] == 0, {i,
nn}], Table[-D[u[xcol[[i]], t], t] +
u[xcol[[i]], t]^kp v2[xcol[[i]], t] +
kp u[xcol[[i]], t]^(kp - 1) u1[xcol[[i]], t] v1[xcol[[i]], t] -
f[xcol[[i]], t] == 0, {i, nn}]];
ic = Join[Table[u[tcol[[i]], 0] == uc[tcol[[i]], 0], {i, nn}],
Table[v[tcol[[i]], 0] == vc[tcol[[i]], 0], {i, nn}]];
bc = {v[x0, t] == vc[x0, t], v[x1, t] == vc[x1, t],
u[x0, t] == uc[x0, t]}; bct = {D[v[x0, t], t] == D[vc[x0, t], t],
D[v[x1, t], t] == D[vc[x1, t], t],
D[u[x0, t], t] == D[uc[x0, t], t]};
var = Join[{u01, v0, v01}, Table[v[i], {i, nn}],
Table[u[i], {i, nn}]]; ic0 = {u[x0, 0] == uc[x0, 0],
v[x0, 0] == vc[x0, 0], v[x1, 0] == vc[x1, 0]};
Finally we solve system of equations and visualize fractional Laplacian (exact and numeric solution) and error in collocation points
sol = NDSolve[{eq, ic, ic0, bct}, var, {t, 0, 1},
Method -> {"EquationSimplification" -> "Residual"}];
{Plot3D[vc[x, t], {x, x0, x1}, {t, 0, 1}, ColorFunction -> "Rainbow",
Mesh -> None, PlotTheme -> "Marketing", AxesLabel -> Automatic,
PlotRange -> All],
Plot3D[v[x, t] /. sol[[1]], {x, x0, x1}, {t, 0, 1},
ColorFunction -> "Rainbow", Mesh -> None, PlotTheme -> "Marketing",
AxesLabel -> Automatic, PlotRange -> All],
Plot[Evaluate[Table[vc[x, t] - v[x, t] /. sol, {x, tcol}]], {t, 0,
1}, PlotLegends -> tcol, AxesLabel -> Automatic, PlotRange -> All]}

It looks nice, but we need to compute $\partial_x v$, and it is why error increases drastically when we using $\partial_x v$ to compute u
lst = Table[{{x, t}, v1[x, t] /. sol[[1]]}, {x,
Join[{x0}, xcol, {x1}]}, {t, 0, 1, .01}];
vx = Interpolation[Flatten[lst, 1], InterpolationOrder -> 4];
U = NDSolveValue[{-D[w[x, t], t] + D[w[x, t]^kp vx[x, t], x] ==
f[x, t], w[x, 0] == uc[x, 0], w[x0, t] == uc[x0, t],
w[x1, t] == uc[x1, t]}, w, {x, x0, x1}, {t, 0, 1}];
Visualization exact (left), numerical solution (center), and error on collocation points
{Plot3D[uc[x, t], {x, x0, x1}, {t, 0, 1}, ColorFunction -> "Rainbow",
Mesh -> None, PlotTheme -> "Marketing", AxesLabel -> Automatic,
PlotRange -> All],
Plot3D[U[x, t], {x, x0, x1}, {t, 0, 1}, ColorFunction -> "Rainbow",
Mesh -> None, PlotTheme -> "Marketing", AxesLabel -> Automatic,
PlotRange -> All],
Plot[Evaluate[Table[uc[x, t] - U[x, t], {x, tcol}]], {t, 0, 1},
PlotLegends -> tcol, AxesLabel -> Automatic, PlotRange -> All]}
We also can ignore last step and compute u
with using sol
. In this case it looks even better then U
{Plot3D[uc[x, t], {x, x0, x1}, {t, 0, 1}, ColorFunction -> "Rainbow",
Mesh -> None, PlotTheme -> "Marketing", AxesLabel -> Automatic,
PlotRange -> All],
Plot3D[u[x, t] /. sol[[1]], {x, x0, x1}, {t, 0, 1},
ColorFunction -> "Rainbow", Mesh -> None, PlotTheme -> "Marketing",
AxesLabel -> Automatic, PlotRange -> All],
Plot[Evaluate[Table[uc[x, t] - u[x, t] /. sol, {x, tcol}]], {t, 0,
1}, PlotLegends -> tcol, AxesLabel -> Automatic, PlotRange -> All]}

Update 1. Note, that in general case the thin film model describes waves. To reproduce wave solutions we can use as initial data ue
and ve
from the paper Finite difference methods for fractional Laplacians in the form
Clear["Global`*"]
ue[x_, k_, s_] := (1 - x^2)^(k + s/2);
lap[x_, k_, s_] :=
Hypergeometric2F1[(1 + s)/2, -k, 1/2, x^2] 2^s Gamma[
k + 1 + s/2] Gamma[(1 + s)/2]/(k! Gamma[1/2]);
The corresponding f
can be evaluated as D[ue[x, k, s]^n D[lap[x, k, s], x], x]
, therefore, we have
f[x_, k_, s_, n_] := (1/(Sqrt[\[Pi]] k!))
2^(2 + s)
k n (k + s/2) (1 + s) x^2 (1 - x^2)^(-1 + k + s/2) ((1 - x^2)^(
k + s/2))^(-1 + n)
Gamma[1 + k + s/2] Gamma[(1 + s)/2] Hypergeometric2F1[1 - k,
1 + (1 + s)/2, 3/2, x^2] - (
2^(1 + s) k (1 + s) ((1 - x^2)^(k + s/2))^
n Gamma[1 + k + s/2] Gamma[(1 + s)/2] Hypergeometric2F1[1 - k,
1 + (1 + s)/2, 3/2, x^2])/(Sqrt[\[Pi]] k!) - (
2^(3 + s) (1 - k) k (1 + s) (1 + (1 + s)/2) x^2 ((1 - x^2)^(
k + s/2))^
n Gamma[1 + k + s/2] Gamma[(1 + s)/2] Hypergeometric2F1[2 - k,
2 + (1 + s)/2, 5/2, x^2])/(3 Sqrt[\[Pi]] k!);
With this function we solve equation $\partial_tu=\partial_x(u^k\partial_x v)-f$ with using code
c[n_, a_] := a 2^(a - 1) Gamma[(a + n)/2]/(Pi^(n/2) Gamma[1 - a/2]);
res = {0 < a < 2};
lapc[n_, a_, x_, u_] :=
c[n, a] Integrate[(u[x] - u[y])/Abs[x - y]^(n + a), y];
UE[m_, t_] := EulerE[m, t]
psi[k_, n_, m_, t_] :=
Piecewise[{{2^(k/2) Sqrt[2/Pi] UE[m, 2^k t - 2 n + 1], (n - 1)/
2^(k - 1) <= t < n/2^(k - 1)}, {0, True}}]
PsiE[k_, M_, t_] :=
Flatten[Table[psi[k, n, m, t], {n, 1, 2^(k - 1)}, {m, 0, M - 1}]]
k0 = 2; M0 = 4; x0 = -1; x1 = 1; nn =
Total[With[{k = k0, M = M0},
Flatten[Table[1, {n, 1, 2^(k - 1)}, {m, 0, M - 1}]]]]; dx = (x1 -
x0)/(nn); xl = Table[x0 + l*dx, {l, 0, nn}]; xcol =
tcol = Table[(xl[[l - 1]] + xl[[l]])/2, {l, 2, nn + 1}]; Psijk =
With[{k = k0, M = M0}, PsiE[k, M, t1]]; Intx1 =
With[{k = k0, M = M0}, Integrate[PsiE[k, M, t1], t1]];
Intx2 = Integrate[Intx1, t1];
intx1[y_] := Intx1 /. t1 -> (y - x0)/(x1 - x0);
intx2[y_] := Intx2 /. t1 -> (y - x0)/(x1 - x0);
Psi[y_] := Psijk /. t1 -> (y - x0)/(x1 - x0);
varu[t_] := Table[u[i][t], {i, nn}];
varv[t_] := Table[v[i][t], {i, nn}];
v[x_, t_] := varv[t] . intx2[x] + x v0[t] + v01[t];
v1[x_, t_] := varv[t] . intx1[x] + v0[t];
v2[x_, t_] := varv[t] . Psi[x];
u[x_, t_] := varu[t] . intx1[x] + u01[t];
u1[x_, t_] := varu[t] . Psi[x];
s = 1/8; kp = 2; L = 2 (1 + s) c[1, s]; mL =
Gamma[2 + s/2] Gamma[5/2 + s/2]; uc[x_] := 1.5 ue[x, 1, s];
vc[x_] := lap[x, 1, s];
int = Table[
Table[NIntegrate[(intx1[tcol[[i]]][[j]] -
intx1[y][[j]])/(tcol[[i]] - y)^(1 + s), {y, x0, tcol[[i]]},
Method -> "PrincipalValue", Exclusions -> tcol[[i]] - y == 0,
AccuracyGoal -> 6, PrecisionGoal -> 6], {j, nn}], {i, nn}] //
Quiet;
int1 = Table[
Table[NIntegrate[(intx1[tcol[[i]]][[j]] -
intx1[y][[j]])/(y - tcol[[i]])^(1 + s), {y, tcol[[i]], x1},
Method -> "PrincipalValue", Exclusions -> y - tcol[[i]] == 0,
AccuracyGoal -> 6, PrecisionGoal -> 6], {j, nn}], {i, nn}] //
Quiet;
int0 = Table[
Table[NIntegrate[(intx1[tcol[[i]]][[
j]])/(tcol[[i]] - y)^(1 + s), {y, -Infinity, x0},
AccuracyGoal -> 10, PrecisionGoal -> 8], {j, nn}], {i, nn}] //
Quiet;
int2 = Table[
Table[NIntegrate[(intx1[tcol[[i]]][[
j]])/(y - tcol[[i]])^(1 + s), {y, x1, Infinity},
AccuracyGoal -> 10, PrecisionGoal -> 8], {j, nn}], {i, nn}] //
Quiet;
(*intb0=Table[NIntegrate[uc[y,0]/(tcol[[i]]-y)^(1+s),{y,-Infinity,x0}]\
,{i,nn}]//Quiet;
intb1=Table[NIntegrate[uc[y,0]/(y-tcol[[i]])^(1+s),{y,x1,Infinity}],{\
i,nn}];*)
lp = c[1, s] (int + int1 + int0 + int2) // Re;
(*eq1=v[x,t]==(-\[CapitalDelta])^(s/2)u[x,t];
eq2=-D[u[x,t],t]+u[x,t]^k D[v[x,t],x,x]+k u[x,t]^(k-1) \
D[u[x,t],x]D[v[x,t],x]+L x D[u[x,t],x]+L u[x,t];*)
eq = Join[
Table[D[varu[t], t] . lp[[i]] - D[v[xcol[[i]], t], t] == 0, {i,
nn}], Table[-D[u[xcol[[i]], t], t] +
u[xcol[[i]], t]^kp v2[xcol[[i]], t] +
kp u[xcol[[i]], t]^(kp - 1) u1[xcol[[i]], t] v1[xcol[[i]], t] -
f[xcol[[i]], 1, s, kp] == 0, {i, nn}]];
ic = Join[Table[u[tcol[[i]], 0] == uc[tcol[[i]]], {i, nn}],
Table[v[tcol[[i]], 0] == vc[tcol[[i]]], {i, nn}]];
bc = {v[x0, t] == vc[x0, t], v[x1, t] == vc[x1, t],
u[x0, t] == uc[x0, t]}; bct = {D[v[x0, t], t] == 0,
D[v[x1, t], t] == 0, D[u[x0, t], t] == 0};
var = Join[{u01, v0, v01}, Table[v[i], {i, nn}],
Table[u[i], {i, nn}]]; ic0 = {u[x0, 0] == 0, v[x0, 0] == vc[x0],
v[x1, 0] == vc[x1]};
sol = NDSolve[{eq, ic, ic0, bct}, var, {t, 0, 1},
Method -> {"EquationSimplification" -> "Residual"}];
Visualization numerical solution with very clear waves for u
(left) and v
(right)
{Plot3D[u[x, t] /. sol[[1]], {x, x0, x1}, {t, 0, 1},
ColorFunction -> "Rainbow", Mesh -> None, PlotTheme -> "Marketing",
AxesLabel -> Automatic, PlotRange -> All, PlotPoints -> 50],
Plot3D[v[x, t] /. sol[[1]], {x, x0, x1}, {t, 0, 1},
ColorFunction -> "Rainbow", Mesh -> None, PlotTheme -> "Marketing",
AxesLabel -> Automatic, PlotRange -> All, PlotPoints -> 50]}

Update 2. We can exclude step with fractional Laplacian computation and compute numerical solution in the first case as follows
c[n_, a_] := a 2^(a - 1) Gamma[(a + n)/2]/(Pi^(n/2) Gamma[1 - a/2]);
lap[n_, a_, x_, u_] :=
c[n, a] Integrate[(u[x] - u[y])/Abs[x - y]^(n + a), y];
ue[a_, x_] := (1 + x^2)^(-(1 - a)/2);
uc[x_, t_] := ue[s, x] (c0 kp (1/(c0 kp) - t))^(-1/kp);
vc[x_, t_] := (c0 kp (1/(c0 kp) - t))^(-1/kp)
2^s Gamma[(1 + s)/2]/Gamma[(1 - s)/2] (1 + x^2)^(-(1 + s)/2);
lape[a_, x_, t_] := (c0 kp (1/(c0 kp) - t))^(-1/kp)
2^a Gamma[(1 + a)/2]/Gamma[(1 - a)/2] (1 + x^2)^(-(1 + a)/2);
f[x_, t_] :=
1/Gamma[(1 - s)/2] (1 - c0 kp t)^(-1 - 1/kp) (1 + x^2)^(-(5/2) - s/
2) (-c0 (1 + x^2)^(2 + s) Gamma[(1 - s)/2] +
2^s kp (-1 + s) (1 + s) (-1 +
c0 kp t) x^2 ((1 - c0 kp t)^(-1/kp) (1 + x^2)^(1/2 (-1 + s)))^
kp Gamma[(1 + s)/2] -
2^s (1 + s) (3 + s) (-1 +
c0 kp t) x^2 ((1 - c0 kp t)^(-1/kp) (1 + x^2)^(1/2 (-1 + s)))^
kp Gamma[(1 + s)/2] +
2^s (1 + s) (-1 + c0 kp t) (1 +
x^2) ((1 - c0 kp t)^(-1/kp) (1 + x^2)^(1/2 (-1 + s)))^
kp Gamma[(1 + s)/2]);
UE[m_, t_] := EulerE[m, t]
psi[k_, n_, m_, t_] :=
Piecewise[{{2^(k/2) Sqrt[2/Pi] UE[m, 2^k t - 2 n + 1], (n - 1)/
2^(k - 1) <= t < n/2^(k - 1)}, {0, True}}]
PsiE[k_, M_, t_] :=
Flatten[Table[psi[k, n, m, t], {n, 1, 2^(k - 1)}, {m, 0, M - 1}]]
k0 = 2; M0 = 4; x0 = 0; x1 = 1; nn =
Total[With[{k = k0, M = M0},
Flatten[Table[1, {n, 1, 2^(k - 1)}, {m, 0, M - 1}]]]]; dx = (x1 -
x0)/(nn); xl = Table[x0 + l*dx, {l, 0, nn}]; xcol =
tcol = Table[(xl[[l - 1]] + xl[[l]])/2, {l, 2, nn + 1}]; Psijk =
With[{k = k0, M = M0}, PsiE[k, M, t1]]; Intx1 =
With[{k = k0, M = M0}, Integrate[PsiE[k, M, t1], t1]];
Intx2 = Integrate[Intx1, t1];
intx1[y_] := Intx1 /. t1 -> (y - x0)/(x1 - x0);
intx2[y_] := Intx2 /. t1 -> (y - x0)/(x1 - x0);
Psi[y_] := Psijk /. t1 -> (y - x0)/(x1 - x0);
varu[t_] := Table[u[i][t], {i, nn}];
u[x_, t_] := varu[t] . intx2[x] + u01[t] + u0[t] x;
u1[x_, t_] := varu[t] . intx1[x] + u0[t];
u2[x_, t_] := varu[t] . Psi[x];
s = 9/10; kp = 3; c0 = 1/6;
intu1 = Table[
Table[NIntegrate[(intx1[tcol[[i]]][[j]] -
intx1[y][[j]])/(tcol[[i]] - y)^(1 + s), {y, x0, tcol[[i]]},
Method -> "PrincipalValue", Exclusions -> tcol[[i]] - y == 0,
AccuracyGoal -> 6, PrecisionGoal -> 6], {j, nn}], {i, nn}] //
Quiet;
int1u1 = Table[
Table[NIntegrate[(intx1[tcol[[i]]][[j]] -
intx1[y][[j]])/(y - tcol[[i]])^(1 + s), {y, tcol[[i]], x1},
Method -> "PrincipalValue", Exclusions -> y - tcol[[i]] == 0,
AccuracyGoal -> 6, PrecisionGoal -> 6], {j, nn}], {i, nn}] //
Quiet;
int0u1 = Table[
Table[NIntegrate[(intx1[tcol[[i]]][[
j]])/(tcol[[i]] - y)^(1 + s), {y, -Infinity, x0},
AccuracyGoal -> 10, PrecisionGoal -> 8], {j, nn}], {i, nn}] //
Quiet;
int2u1 = Table[
Table[NIntegrate[(intx1[tcol[[i]]][[
j]])/(y - tcol[[i]])^(1 + s), {y, x1, Infinity},
AccuracyGoal -> 10, PrecisionGoal -> 8], {j, nn}], {i, nn}] //
Quiet;
intu2 = Table[
Table[NIntegrate[(Psi[tcol[[i]]][[j]] -
Psi[y][[j]])/(tcol[[i]] - y)^(1 + s), {y, x0, tcol[[i]]},
Method -> "PrincipalValue", Exclusions -> tcol[[i]] - y == 0,
AccuracyGoal -> 6, PrecisionGoal -> 6], {j, nn}], {i, nn}] //
Quiet;
int1u2 = Table[
Table[NIntegrate[(Psi[tcol[[i]]][[j]] -
Psi[y][[j]])/(y - tcol[[i]])^(1 + s), {y, tcol[[i]], x1},
Method -> "PrincipalValue", Exclusions -> y - tcol[[i]] == 0,
AccuracyGoal -> 6, PrecisionGoal -> 6], {j, nn}], {i, nn}] //
Quiet;
int0u2 = Table[
Table[NIntegrate[(Psi[tcol[[i]]][[
j]])/(tcol[[i]] - y)^(1 + s), {y, -Infinity, x0},
AccuracyGoal -> 10, PrecisionGoal -> 8], {j, nn}], {i, nn}] //
Quiet;
int2u2 = Table[
Table[NIntegrate[(Psi[tcol[[i]]][[
j]])/(y - tcol[[i]])^(1 + s), {y, x1, Infinity},
AccuracyGoal -> 10, PrecisionGoal -> 8], {j, nn}], {i, nn}] //
Quiet;
lpu1 = c[1, s] (intu1 + int1u1 + int0u1 + int2u1) // Re;(*v1*)
lpu2 = c[1, s] (intu2 + int1u2 + int0u2 + int2u2) // Re;(*v2*)
eq = Table[-D[u[xcol[[i]], t], t] +
u[xcol[[i]], t]^kp varu[t] . lpu2[[i]] +
kp u[xcol[[i]], t]^(kp - 1) u1[xcol[[i]], t] varu[t] .
lpu1[[i]] - f[xcol[[i]], t] == 0, {i, nn}];
ic = Table[u[tcol[[i]], 0] == uc[tcol[[i]], 0], {i, nn}];
bc = {u[x0, t] == uc[x0, t],
u[x1, t] == uc[x1, t]}; bct = {D[u[x1, t], t] == D[uc[x1, t], t],
D[u[x0, t], t] == D[uc[x0, t], t]};
var = Join[{u01, u0},
Table[u[i], {i, nn}]]; ic0 = {u[x0, 0] == uc[x0, 0],
u[x1, 0] == uc[x1, 0]};
sol = NDSolve[{eq, ic, ic0, bct}, var, {t, 0, 1},
Method -> {"EquationSimplification" -> "Residual"}];
Visualization exact and numerical solution, and error in the collocation points
{Plot3D[uc[x, t], {x, x0, x1}, {t, 0, 1}, ColorFunction -> "Rainbow",
Mesh -> None, PlotTheme -> "Marketing", AxesLabel -> Automatic,
PlotRange -> All],
Plot3D[u[x, t] /. sol[[1]], {x, x0, x1}, {t, 0, 1},
ColorFunction -> "Rainbow", Mesh -> None, PlotTheme -> "Marketing",
AxesLabel -> Automatic, PlotRange -> All],
Plot[Evaluate[Table[uc[x, t] - u[x, t] /. sol, {x, tcol}]], {t, 0,
1}, PlotLegends -> tcol, AxesLabel -> Automatic, PlotRange -> All]}
