# differential equations – Why my code give empty plots? AND When I increase the number of mx , mt my code is running for a long time

In my code, two graphs are not displayed. what is the problem?Also when I increase the number of mx and mt, the code enters the running for a long time.
The regularization problem is described as
$$D^{alpha}w_t(x,t)=w_{xx}(x,t) , 0 , (1)

$$w(x,0)=J_{delta_{1}}phi^epsilon(x),0leq xleq s(0),$$(2)
$$w(s(t),t)=0,0leq tleq T,$$(3)

$$w_(s(t),t)=-J_{delta_2}g^epsilon_2(t),0leq tleq T,$$(4)

where $$J_delta(.)$$shows the mollification function with respect to the radius of mollification,$$delta$$.

``````(*Inverse fractional Stefan problem by mollification and Marching
scheme*)
ClearAll(a, exact, g, g1, g2, f, epsilon, deltax);
a = .5;(*a= alpha*)
T = 1;
exact(x_, t_) = .5*Sqrt(Pi)*x^2 + 2 Sqrt(t);(*   u(x,t)  *)
phi(x_) = exact(x, 0);(*   u(x,0)=(Phi)  *)
g(t_) = Sqrt(t^2 + 2);(*   s(t)  *)
g1(t_) = exact(g(t), t);(*   u(s(t),t)=Subscript(g, 1)(t)  *)
g2(t_) = D(exact(x, t), x) /.
x -> g(t);(*   Subscript(u, x)(s(t),t)=Subscript(g, 2)(t)  *)
f(s_, z_) = Sqrt(Pi)*.5 - D(exact(s, z), {s, 2});
epsilon = 0.01*Sqrt(3);
deltax = 0.05; deltat = 0.07;
Mx = 20;
Mt = 30;
h = N(g(0)/Mx);
k = N(T/Mt);
rx = Floor((3*deltax)/h);
rt = Floor((3*deltat)/k);
Array(tim, Mt + 3, -1);
Array((Theta), Mt + 1, 0);
Array(M, Mt + 1, 0);
Do(tim(j) = j*k, {j, -1, Mt + 1});
Do(M(n) = Floor(g(tim(n))/h); (Theta)(n) = g(tim(n))/h - M(n), {n, 0,
Mt});
Do(points(i) = Array(meshpoint, {1, M(i) + 1}, {i, 0}), {i, 0, Mt});
Do(points1(i) = Array(meshpoint1, {1, M(i) + 1}, {i, 0}), {i, 0, Mt});
Do(meshpoint(i, j) = j*h, {i, 0, Mt}, {j, 0, M(i)})
Array(points, Mt + 1, 0);
MatrixForm(Array(points, Mt + 1, 0));
Do(meshpoint1(i, j) = {j*h, tim(i)}, {i, 0, Mt}, {j, 0, M(i)});
MatrixForm(Array(points1, Mt + 1, 0));
ListPlot(Flatten(Array(points1, Mt + 1, 0), 2));
Ap = 1/NIntegrate(Exp(-t^2), {t, -3, 3});
ro(t_, delta_) = Ap/delta*Exp(-t^2/delta^2);
Do(s2(i) = i*h, {i, -rx, Mx + rx});
Do(s2t(i) = i*k, {i, -rt, Mt + rt});
a12 = epsilon*RandomReal({-1, 1}, {Mx + 1});
Do(noisephi(i) = phi(meshpoint(0, i)) + a12((i + 1)), {i, 0, Mx});
RO(x_) = g(0)/(Mx + 1) Sum(
ro(x - s2(i), deltax)*noisephi(i), {i, 0, Mx});
c1 = NIntegrate((phi(x) - (RO(x)))*
Integrate(ro(x - z, deltax), {z, -3*deltax, 0}), {x, 0,
3*deltax})/
NIntegrate(
Integrate(ro(x - z, deltax), {z, -3*deltax, 0})^2, {x, 0,
3*deltax});
c2 = NIntegrate((phi(x) - (RO(x)))*
Integrate(ro(x - z, deltax), {z, g(0), g(0) + 3*deltax}), {x,
g(0) - 3*deltax, g(0)})/
NIntegrate(
Integrate(ro(x - z, deltax), {z, g(0), g(0) + 3*deltax})^2, {x,
g(0) - 3*deltax, g(0)});
If(rx > 0,
RO1(x_) = (
3*deltax)/(rx + 1) Sum(ro(x - s2(i), deltax)*c1, {i, -rx, -1}),
RO1(x_) = 0);
RO3(x_) = (
3*deltax)/(rx + 1) Sum(
ro(x - s2(i), deltax)*c2, {i, Mx + 1, Mx + rx});
jphi2(x_) =
Piecewise({{RO1(x) + RO(x), 0 <= x < 3*deltax}, {RO(x),
3*deltax <= x < g(0) - 3*deltax}, {RO3(x) +
RO(x), (g(0) - 3*deltax) <= x <= g(0)}});
(1/(Mx + 1)
Sum((exact(meshpoint(0, j), 0) - jphi2(meshpoint(0, j)))^2, {j, 0,
Mx}))^(1/2)
o1 = Plot(phi(x), {x, 0, g(0)})
o2 = Plot(jphi2(x), {x, 0, g(0)})
Show(o1, o2)
(**********************jg2*****************************************)
Do(
noiseg2(n) = g2(n*k) + a12((n + 1)), {n, 0, Mt});
ROg2(t_) =
1/(Mt + 1) Sum(ro(t - s2t(i), deltat)*noiseg2(i), {i, 0, Mt});
c1g2 = NIntegrate((g2(t) - (ROg2(t)))*
Integrate(ro(t - z, deltat), {z, -3*deltat, 0}), {t, 0,
3*deltat})/
NIntegrate(
Integrate(ro(t - z, deltat), {z, -3*deltat, 0})^2, {t, 0,
3*deltat});
c2g2 = NIntegrate((g2(t) - (ROg2(t)))*
Integrate(ro(t - z, deltat), {z, 1, 1 + 3*deltat}), {t,
1 - 3*deltat, 1})/
NIntegrate(
Integrate(ro(t - z, deltat), {z, 1, 1 + 3*deltat})^2, {t,
1 - 3*deltat, 1});
If(rt > 0,
RO1g2(t_) = (
3*deltat)/(rt + 1) Sum(ro(t - s2t(i), deltat)*c1g2, {i, -rt, 0}),
RO1(t_) = 0);
RO3g2(t_) = (
3*deltat)/(rt + 1) Sum(
ro(t - s2t(i), deltat)*c2g2, {i, Mt, Mt + rt});
RO4g2(t_) = (
1 + 3*deltat)/(Mt + rt + 1) Sum(
ro(t - s2t(i), deltat)*noiseg2(i), {i, 0, Mt});
jg2(t_) =
Piecewise({{RO1g2(t) + ROg2(t), 0 <= t < 3*deltat}, {ROg2(t),
3*deltat <= t < 1 - 3*deltat}, {RO3g2(t) +
RO4g2(t), (1 - 3*deltat) <= t <= 1}});
o1g = Plot(g2(t), {t, 0, 1})
o2g = Plot(jg2(t), {t, 0, 1})
Show(o1g, o2g)
(1/(Mt + 1) Sum(((g2(n*k)) - jg2(n*k))^2, {n, 0, Mt}))^(1/2)

gt(x_) = If(x >= g(0), (3 - (2 - x)^2)/2, 0);
Array(num, M(Mt) + 1, 0);
Do(num(i) =
If(gt(i*h) < 0, 0,
If(gt(i*h)/k == Floor(gt(i*h)/k), gt(i*h)/k,
Floor(gt(i*h)/k) + 1)), {i, 0, M(Mt)});
Array(num, M(Mt) + 1, 0)
(*************************************************************************************************)

Do(
u1(j, 0) = jphi2(j*h), {j, 0, M(0)});
Q1(0, 0) = (jphi2(h) - jphi2(0))/(h);
Q1(M(0), 0) = (jphi2((M(0))*h) - jphi2((M(0) - 1)*h))/(h);
Do(
Q1(j, 0) = (jphi2((j + 1)*h) - jphi2((j - 1)*h))/(2*h), {j, 1,
M(0) - 1});
Do(
Q1(M(n) + (Theta)(n), n) = jg2(n*k);
u1(M(n) + (Theta)(n), n) = 0;
w1(M(n) + (Theta)(n), n) = -1
, {n, 0, Mt});
Do(
u1(M(n),
n) = -(Theta)(n)*h*Q1(M(n) + (Theta)(n), n) +
u1(M(n) + (Theta)(n), n) // N;
Q1(M(n), n) = -(Theta)(n)*h*(-f((M(n) + (Theta)(n))*h, n*k)) +
Q1(M(n) + (Theta)(n), n);
w1(M(n), n) = w1(M(n) + (Theta)(n), n)
, {n, 1, Mt});
Do(Do(ww(i_) :=
N(1/Gamma(1 - a)*1/(
1 - a)*(((2*i + 1)*deltat/2)^(1 - a) - ((2*i - 1)*deltat/2)^(1 -
a))), {i, 1, j - 1}), {j, 1, Mt});
Do(wwj(jj_) :=
N(1/Gamma(1 - a)*1/(
1 - a)*(jj*deltat - ((jj - .5)*deltat)^(1 - a))), {jj, 1, Mt});
ww(1) = N(1/Gamma(1 - a)*1/(1 - a)*(deltat/2)^(1 - a));

Do(rq3(M(Mt) - j) = Floor(((3*(1 - num(M(Mt) - j)*k))/6)/k);
u1(M(Mt) - j - 1, n) = u1(M(Mt) - j, n) - h*Q1(M(Mt) - j, n) // N;
Q1(M(Mt) - j - 1,
n) = -h*(w1(M(Mt) - j, n) - f((M(Mt) - j)*h, n*k)) +
Q1(M(Mt) - j, n);
jQ1(M(Mt) - j, t_) =
If(num(-j + M(Mt)) == Mt,
0, (1 - num(M(Mt) - j)*k -
6*((1 - num(M(Mt) - j)*k)/6))/(Mt - num(M(Mt) - j) + 1 -
2*rq3(M(Mt) - j)) Sum(
ro(t - s1(d), (1 - num(M(Mt) - j)*k)/6)*Q1(M(Mt) - j, d), {d,
num(M(Mt) - j) + rq3(M(Mt) - j), Mt - rq3(M(Mt) - j)}));

df(M(Mt) - j, i_) :=
N((jQ1(M(Mt) - j, tim(i + 1)) - jQ1(M(Mt) - j, tim(i)))/
deltat);(*forward difference*)
db(M(Mt) - j, i_) :=
N((jQ1(M(Mt) - j, tim(i)) - jQ1(M(Mt) - j, tim(i - 1)))/
deltat);(*backward difference*)
d0(M(Mt) - j, i_) :=
N((jQ1(M(Mt) - j, tim(i + 1)) - jQ1(M(Mt) - j, tim(i - 1)))/(
2*deltat));(*central difference*)
da(M(Mt) - j, 1) =
df(M(Mt) - j, 1)*ww(1)    ;                            (*
Subscript((D^(Alpha)G^(Epsilon)), (Delta))(Subscript(t, 1)) *)
da(M(Mt) - j, 2) =
df(M(Mt) - j, 1)*ww(1) + df(M(Mt) - j, 2)*wwj(2); (*
Subscript((D^(Alpha)G^(Epsilon)), (Delta))(Subscript(t, 2)) *)
Do(da(M(Mt) - j, k) =
df(M(Mt) - j, 1)*ww(1) + df(M(Mt) - j, k)*wwj(k) +
Sum(d0(M(Mt) - j, i)*ww(k - i + 1), {i, 2, k - 1}), {k, 3,
n});(*   Subscript((D^(Alpha)G^(Epsilon)), (Delta))(Subscript(
t, j)) *)

w1(M(Mt) - j - 1, n) = w1(M(Mt) - j, n) - h*da(M(Mt) - j, n)
, {j, 0, M(Mt) - M(0) - 1}, {n, num(-j + M(Mt)), Mt});
Do(
rq = Floor(3*.16666/k);
u1(M(Mt) - j - 1, n) = u1(M(Mt) - j, n) - h*Q1(M(Mt) - j, n) // N;
Q1(M(Mt) - j - 1,
n) = -h*(w1(M(Mt) - j, n) - f((M(Mt) - j)*h, n*k)) +
Q1(M(Mt) - j, n);
jQ1(M(Mt) - j, t_) =
If(num(-j + M(Mt)) == Mt,
0, (1 - 6*.16666)/(Mt - 2*rq + 1) Sum(
ro(t - s1(d), .16)*Q1(M(Mt) - j, d), {d, rq,
Mt - rq}));(*Runing*)
w1(M(Mt) - j - 1, n) = w1(M(Mt) - j, n) - h*da(M(Mt) - j, n + 1)
, {j, M(Mt) - M(0), M(Mt) - 1}, {n, 0, Mt});
fi = Table({tim(i), u1(0, i)}, {i, 0, Mt});
fi1 = ListLinePlot(fi, PlotStyle -> {Dashed, Red}, Frame -> True)
fi2 = Plot(exact(0, t), {t, 0, 1}, PlotStyle -> Blue, Frame -> True)
Show(fi1, fi2, GridLines -> Automatic,
GridLinesStyle -> Directive(Blue, Dotted))
ss3 = Table({tim(n), Abs(u1(0, n) - exact(0, n*k))}, {n, 0, Mt});
gg5 = ListPlot(ss3, PlotStyle -> Red, Frame -> True,
GridLines -> Automatic, GridLinesStyle -> Directive(Blue, Dotted))
``````