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

This notebook focuses on producing the reference values for c(0,t) and (Subscript[x, 0])(t) for t = 0.01, 0.02, ..., 0.19, 0.1974. Calculations are done with stepsizes dx = 0.002, 0.001, 0.0005, 0.00025, 0.000125 and then Richardson extrapolation is used.

prec=30;
Block[{$MinPrecision=prec},
tf=1974/10000; 
tf1=1/100;

wa=1/2;
wb=3/2;
lst={};

For[nx=500,nx<8001,nx=2nx,   
dx=N[1/nx,prec];

α=1/5;dt=α/nx^2; α2=2α;β=1-α2;
nts=tf1/dt;
nt=tf/dt;

c=Table[(1-i dx)^2/2,{i,0,nx}];
r=nx;              (* start with c[[r+1]]=0, corresponding to x_r+1 = r dx *)

t=0;x0=1; dt=N[dt,prec];
xl={{t,x0}};       (* xl is a list of boundary points *)


v1=0;v2=0;v0=0; (* for tracking branches used *)
vt=0;vts=0;
While[vt<nt,
vt++;vts++;

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

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;
(*xl=Append[xl,{t,x0}];*)

y=x0-xr;

If[x0<xr+wa dx, r--;v0++;c=Most[c],                                              
y=y-dx;v1++;
c[[r+1]]=y^2(1/2+y x0p/6 + (y x0p)^2/24);
If[x0>xr+wb dx, r++;c=Append[c,0];v2++;
y=y-dx;
c[[r+1]]=y^2(1/2+y x0p/6 + (y x0p)^2/24)]];

If[vts≥nts,vts=0;lst=Append[lst,{dx,t,c[[1]],x0}];Print[{NumberForm[dx,3],NumberForm[t,2],c[[1]],x0,v0,v1,v2}]]
];
lst=Append[lst,{dx,t,c[[1]],x0}];Print[{NumberForm[dx,3],NumberForm[t,4],c[[1]],x0,v0,v1,v2}]
]];
TableForm[lst]

r
c[[r+1]]

x=0;lstc={};
For[i=1,i<=r+1,i++,
AppendTo[lstc,{x,c[[i]]}];x=x+dx];
ListPlot[lstc]
lstc

ListPlot[lstc,Joined->True,AxesLabel->{"x","c(x,0.1974)"}]

Getting both Richardson extrapolations and error estimates by processing lst

prec=30;
sA[t_]:=2N[Erf[1/(2Sqrt[t])],prec]-1; (* asymptotic formula for s(t) *)
cA[t_]:=N[1/2-2Sqrt[t/Pi],prec];   (* asymptotic formula for c(0,t) *)

steps=20; 
For[i=1,i≤steps,i++, 
tt=lst[[i,2]];    (* time *)
co=lst[[i,3]];so=lst[[i,4]];     (* h *)
lst1={{"dx","  c(0)  ","Rc","ERc/dx^2" ,"R2c","ER2c/dx^4","s","Rs","ERs/dx^2" ,"R2s","ER2s/h^4"},{lst[[i,1]],co,"","","","",so,"","","",""}};

j=i+steps;k=i;         (* h/2 *)
h=lst[[j,1]]; cn=lst[[j,3]];sn=lst[[j,4]];Rc1=(4cn-co)/3;ERc=(cn-co)/(3h^2); 
Rs1=(4sn-so)/3; ERs=(sn-so)/(3h^2);co=cn;so=sn;
AppendTo[lst1,{h,cn,Rc1,ERc,"","",sn,Rs1,ERs,"",""}];      (* 3rd row of lst1 *)

k=j;j=j+steps;     (* h/4 *)
h=lst[[j,1]]; cn=lst[[j,3]];sn=lst[[j,4]];Rc2=(4cn-co)/3;ERc=(cn-co)/(3h^2); 
Rs2=(4sn-so)/3; ERs=(sn-so)/(3h^2);co=cn;so=sn;
R2c=(16Rc2-Rc1)/15; ER2c=(Rc2-Rc1)/(15h^4);
R2s=(16Rs2-Rs1)/15;ER2s=(Rs2-Rs1)/(15h^4);
Rc1=Rc2;Rs1=Rs2; 
AppendTo[lst1,{h,cn,Rc2,ERc,R2c,ER2c,sn,Rs2,ERs,R2s,ER2s}];  (* 4th row of lst1 *)

k=j;j=j+steps;       (* h/8 *)
h=lst[[j,1]]; cn=lst[[j,3]];sn=lst[[j,4]];Rc2=(4cn-co)/3;ERc=(cn-co)/(3h^2); 
Rs2=(4sn-so)/3; ERs=(sn-so)/(3h^2);
R2c=(16Rc2-Rc1)/15; ER2c=(Rc2-Rc1)/(15h^4);
R2s=(16Rs2-Rs1)/15;ER2s=(Rs2-Rs1)/(15h^4);
Rc1=Rc2;Rs1=Rs2; co=cn;so=sn;
AppendTo[lst1,{h,cn,Rc2,ERc,R2c,ER2c,sn,Rs2,ERs,R2s,ER2s}];(* 5th row of lst1 *)

k=j;j=j+steps;       (* h/16 *)
h=lst[[j,1]]; cn=lst[[j,3]];sn=lst[[j,4]];Rc2=(4cn-co)/3;ERc=(cn-co)/(3h^2); 
Rs2=(4sn-so)/3; ERs=(sn-so)/(3h^2);
R2c=(16Rc2-Rc1)/15; ER2c=(Rc2-Rc1)/(15h^4);
R2s=(16Rs2-Rs1)/15;ER2s=(Rs2-Rs1)/(15h^4);
AppendTo[lst1,{h,cn,Rc2,ERc,R2c,ER2c,sn,Rs2,ERs,R2s,ER2s}]; (* 6th row of lst1 *)

If[i<6, (* comparing to assymtotic values *)
cna=cA[tt];sna=sA[tt];
AppendTo[lst1,{"Asym", cna,"",( cna-cn)/h^2,"",(cna-Rc2)/h^4,sna,"",(sna-sn)/h^2,"",(sna-Rs2)/h^4}]];


Print["At t = ",tt,"\n",TableForm[lst1]];

mm[v_]:=MinMax[{lst1[[6,v]],lst1[[5,v]],lst1[[4,v]]}];
range[v_]:=mm[v][[2]]-mm[v][[1]];
ERce=h^2 range[4]; dRc=Ceiling[Log10[Rc2/ERce]+.5];
ER2ce=h^4 range[6];dR2c=Ceiling[Log10[R2c/ER2ce]+.5];
ERse=h^2 range[9]; dRs=Ceiling[Log10[Rs2/ERse]+.5];
ER2se=h^4 range[11];dR2s=Ceiling[Log10[R2s/ER2se]+.5];

Print["Rc = ",Rc2,"  from ", Rc2-ERce," to ", Rc2+ERce," ** error = ",ERce];     
Print["Rc = ", NumberForm[Rc2,dRc],"  from ",NumberForm[Rc2-ERce,dRc]," to ", NumberForm[Rc2+ERce,dRc]," ** error = ",NumberForm[ERce,3]];  
Print["R2c = ",R2c,"  from ", R2c-ER2ce," to ", R2c+ER2ce," ** error = ",ER2ce];
Print["R2c = ",NumberForm[R2c,dR2c],"  from ", NumberForm[R2c-ER2ce,dR2c]," to ", NumberForm[R2c+ER2ce,dR2c]," ** error = ",NumberForm[ER2ce,3]];
Print["Rs = ",Rs2,"  from ",Rs2-ERse," to ",Rs2 + ERse, " ** error = ",ERse];
Print["Rs = ",NumberForm[Rs2,dRs],"  from ",NumberForm[Rs2-ERse,dRs]," to ",NumberForm[Rs2 + ERse,dRs], " ** error = ",NumberForm[ERse,3]];
Print["R2s = ",R2s,"  from ",R2s-ER2se," to ",R2s + ER2se," ** error = ",ER2se];
Print["R2s = ",NumberForm[R2s,dR2s],"  from ",NumberForm[R2s-ER2se,dR2s]," to ",NumberForm[R2s + ER2se,dR2s], " ** error = ",NumberForm[ER2se,3],"\n"]
];