Crank&Gupta, A Moving boundary problem
((∂c)/(∂t))=(((∂^2)c)/(∂(x^2)))-1      for   0 < x < (Subscript[x, 0])(t),    t>0
((∂c)/(∂x))(0,t)=0,             for t>0
((∂c)/(∂x))((Subscript[x, 0])(t),t)=c((Subscript[x, 0])(t),t)=0,   for t>0
u(x,0)=(1-x)2/2           for    0<x<1

Continuation from given c,s given in newCG01ffa2Prec30xxxxx (10/13/22),  r=223, t=0.1974
to t=0.197434+

data=lstc

prec=30;
Block[{$MinPrecision=prec},
ts=1974/10000; (* starting time for continuation *)
x0=0.0278290124988106673770944286210370522618186564699736132904`30.;

wa=1/2; wb=3/2;  (* determine the 'no change window' *)



dx=1/8000;
r=223;
α=1/5;α2=2α;β=1-α2;
dt=α dx^2;
nt=(197435/1000000-ts)/dt-8;

c=Table[data[[i,2]],{i,1,r+1}];   (* initial values of the concentration *)
Print[nt, "  ",Dimensions[c]];

t=ts; dt=N[dt,prec];     (* position of the boundary at t is x0 = s(t)*)
xl={{t,x0}};
cl={{t,c[[1]]}};

vt=0;
While[vt< nt,
vt++;

cn=β c;
cn[[1]]+= α2 c[[2]]-dt; (* for 0 flux *)
t+=dt;
AppendTo[cl,{t,c[[1]]}];

For[i=2,i≤r,i++, 
cn[[i]]+= α (c[[i+1]]+c[[i-1]])-dt];
c=cn;

xr=(r-1)dx;

y1=Sqrt[2c[[r]]];
x0p=(y1-x0+xr)/(dt+y1^2/6); 
x0+=x0p dt;
AppendTo[xl,{t,x0}];
y=x0-xr;

If[x0<xr+wa dx, r--;c=Most[c],                                              
y=y-dx;
c[[r+1]]=y^2(1/2+y x0p/6 + (y x0p)^2/24);
If[x0>xr+wb dx, r++;c=Append[c,0]; (* this branch is not used because s is only decreasing for the classical oxygen diffusion problem*)
y=y-dx;
c[[r+1]]=y^2(1/2+y x0p/6 + (y x0p)^2/24)]]
];
Print[t,"  ",r,"  ", c[[r+1]]];
x=0;lstc={};
For[i=1,i<=r+1,i++,
AppendTo[lstc,{x,c[[i]]}];x=x+dx];
Print[lstc]];

ListPlot[lstc]
xl
cl

n=Dimensions[cl][[1]];
last=9;
lc=Table[{cl[[i,2]],cl[[i,1]]},{i,n-last,n}];
lc0=ListPlot[lc,PlotStyle->Red,Joined->False];
t=.;
linec[t_]=Fit[lc,{1,t},t];
tc=linec[0];
Print["Extiction time using last ",last+1," points of c(0,t) is ",NumberForm[tc,10]]
lc2=Plot[linec[t],{t,0,0.00000005},PlotStyle->{Black,Thin},AxesLabel->{"    c(0,t)","t"}];
Show[lc2,lc0]

last=499;
lc=Table[{cl[[i,2]],cl[[i,1]]},{i,n-last,n}];
lc0=ListPlot[lc,PlotStyle->Red,Joined->False];
lc2=Plot[linec[t],{t,0,0.000002},PlotStyle->{Black,Thin},AxesLabel->{"    c(0,t)","t"}];
Show[lc2,lc0]

skip=200; (* skip last 200 points *)
last=9;
lcs=Table[{cl[[i,2]],cl[[i,1]]},{i,n-last-skip,n-skip}];
lcs0=ListPlot[lcs,PlotStyle->Red,Joined->False];
t=.;
linecs[t_]=Fit[lcs,{1,t},t];
tcs=linecs[0];
Print["Skip last ",skip," points. Extiction time using the last ",last+1," of the remaining points of c(0,t) is ",NumberForm[tcs,10]];
lcs=Plot[{linec[t],linecs[t]},{t,0,0.0000008},(*PlotStyle->{Black,Thin},*)AxesLabel->{"    c(0,t)","t"}];
Show[lcs,lcs0]

last=290;skip=20;
lca=Table[{(tc-cl[[i,1]])^.25,cl[[i,2]]/(tc-cl[[i,1]])},{i,n-last-skip,n-skip}];
lc0a=ListPlot[lca,PlotStyle->Red,Joined->False]

n=Dimensions[xl][[1]];
last=9;factor=10000;
end=Table[{factor xl[[i,2]]^2,xl[[i,1]]},{i,n-last,n}];
g2=ListPlot[end,PlotStyle->Red,Joined->False];

t=.;
lines[t_]=Fit[end,{1,t},t];
tcs=lines[0];
Print["Extiction time using last ",last," points of s is ",NumberForm[tcs,10]]

g1=Plot[lines[t],{t,0,.012},PlotStyle->{Black,Thin}, AxesLabel->{"   (100 s(t))^2 ","t"}];
Show[g1,g2]

last=499;factor=10000;
end=Table[{factor xl[[i,2]]^2,xl[[i,1]]},{i,n-last,n}];
g2=ListPlot[end,PlotStyle->Red,Joined->False];
g1=Plot[lines[t],{t,0,0.4},PlotStyle->{Black,Thin}, AxesLabel->{"   (100 s(t))^2 ","t"}];
Show[g1,g2]

skip=200; (* skip last 200 points *)
last=9;
ldss=Table[{factor xl[[i,2]]^2,xl[[i,1]]},{i,n-last-skip,n-skip}];
lss0=ListPlot[ldss,PlotStyle->Red,Joined->False];
t=.;
liness[t_]=Fit[ldss,{1,t},t];
tss=liness[0];
Print["Skip last ",skip," points. Extiction time using the last ",last+1," of the remaining points of s(t) is ",NumberForm[tss,10]];
lss=Plot[{lines[t],liness[t]},{t,0,0.17},(*PlotStyle->{Black,Thin},*)AxesLabel->{"   (100 s(t))^2 ","t"}];
Show[lss,lss0]

last=5999;skip=5;
lsa=Table[{((tcs-xl[[i,1]]))^.1,xl[[i,2]]^2/(-xl[[i,1]]+tcs)},{i,n-last-skip,n-skip}];
ls0a=ListPlot[lsa,PlotStyle->Red,Joined->False];

t=.;
linesx[t_]=Fit[lsa,{1,t},t];
tcxs=linesx[0];
NumberForm[tcxs,10]

g1x=Plot[linesx[t],{t,0.,.4},PlotStyle->{Black,Thin}, AxesLabel->{"(tc-t)^.1"," s(t)^2/(tc-t) "}];
Show[g1x,ls0a]