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. Calculations are done with stepsizes dx = 0.002, 0.001, 0.0005, 0.00025, 0.000125 and then Richardson extrapolation is to be used.

prec=30;
Block[{$MinPrecision=prec},
tf=1/100; (* final time, increase it if you want to, I have done it with 1974/10000 *)
tf1=1/100; (* output at multiples of tf1 *)

wa=1/2; wb=3/2;  (* determine the 'no change window' *)
lst={};                 (* output *)
Print["   dx      t       c(0,t)                           s(t)"];


For[nx=500,nx<8001,nx=2nx,   
dx=N[1/nx,prec];    (* calulations will be done with different dx *)

α=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}];   (* initial values of the concentration *)
r=nx;              (* start with c[[r+1]]=0, corresponding to x_r+1 = r dx *)

t=0;x0=1; dt=N[dt,prec];     (* position of the boundary at t is x0 = s(t)*)

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;

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)]];

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