Classical moving boundary problem
((∂c)/(∂t))=(((∂^2)c)/(∂(x^2)))-1      for   0 < x < s(t),    t>0
((∂c)/(∂x))(0,t)=0            for t>0     (FLUX =0)
((∂c)/(∂x))(s(t),t)=c(s(t),t)=0,   for t>0
c(x,0)=(1-x)^2/2,   for    0<x<1

(* compuatation of the approximations by the numerical method *)
tf=.16;       (* final time *)


wa=0.3;
wb=1.31;
Print["Window [ ",wa," , ",wb," ]dx"];

nx=320;   
Print["dx = ",1/nx];
dx=1./nx;  

α=.4;dt=α dx^2;
nt=Ceiling[tf/dt];dt=tf/nt;α=dt/dx^2;
Print["α = ",α];

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

t=0.;x0=1.;
xl={{t,x0}};       (* xl is a list of boundary points *)

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

cn=(1-2α)c;
cn[[1]]+= 2α( c[[2]](*-dx flux[t]*))-dt;     (* for prescribed 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=(2c[[r]])^.5;
x0p=(y1-x0+xr)/(dt+y1^2/6); (* drop y1^2/6 to get Crank-Gupta placement *)
x0+=x0p dt;
xl=Append[xl,{t,x0}];

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

];           (* reached final t *)




Print["At t = ",tf," approximate s = ",x0];
Print["At x = 0, t = ",tf," approximate c = ",c[[1]]];