A moving boundary problem
((∂c)/(∂t))=(((∂^2)c)/(∂(x^2)))-1      for   0 < x < s(t),    t>0
((∂c)/(∂x))(s(t),t)=c(s(t),t)=0,   for t>0
((∂c)/(∂x))(0,t)=flux(t)            for t>0     (FLUX PRESCRIBED)
c(x,0)=initial(x)           for    0<x<1
In the next cell specify functions flux and initial and run it.
The cell following will compute the solution at final time tf

(* Example 2: for prescribed s(t)=(1-5t)^.5 analytic solution of the moving boundary problem is found in this cell, which determines flux(t) and initial(x). The numerical method, with this flux(t) and initial(x), is used to recover s(t) and c(x,t) in the following cell. *)

t0=1/5;
a=2t0;
w[z_]:=Module[{s=z^2/2,a0=1,a1=-1,b,t=1,f=z/3,n=4},
f*=s; s-=f;
While[Abs[t]>10^-17&&n<100,b=-a1+a (n-4) a0;a0=a1;a1=b;
f*=z /n;t=f b;s+=t;n++];
s];
u[x_,t_]:=a^2(1-t/t0)w[(1-x/(1-t/t0)^.5)/a];     (* analytict solution *)
initial[x_] := u[x,0];

a0=1;a1=-1.;
z1=1./(2a^2); c01=1/a-z1;n=3;
While[Abs[z1 a1]>10^-17 && n<100, 
b=-a1+a (n-3) a0;a0=a1;a1=b;z1=z1/(n a);c01+=z1 b;n++];
Print["Number of terms used = ",n];
c01=-2c01 t0^.5;   (* = -2 w'(1/a) t0^.5*)
flux[t_] := c01 (t0-t)^.5;

(* DIVERSION (=not needed) : Integral computation of c[x,t]: *)
xt=.5;tt=.1;st=(1-tt/t0)^.5;
Print["Analytic value of c(x,t) at x = ",xt,",  t = ",tt,":  ",u[xt,tt]];
Print["Integral version:  ",
st^2/2((1-xt/st)^2+(1/a)NIntegrate[(y-xt/st)^2E^((y^2-1)/(2a)),{y,1,xt/st}])];

(* compuatation of the approximate s(t) and c(x,t) by the numerical method for given flux(t) and initial(x)*)
tf=.16;       (* final time *)

(* at tf to be compared with:*)
x0f=(1-tf/t0)^.5;    (* exact value of s(tf) *)
rx=2/10;                           (* reporting x at the end, enter it as fraction! *)
vf=u[rx,tf];                 (* exact value of c(rx,tf) *)


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[initial[i dx],{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[ListPlot[Table[{i dt, xl[[i+1]][[2]]-(1-i dt/t0)^.5},{i,0,nt}],PlotRange->All,
AxesLabel->{"t","Error s(t)"}]];

cx=c[[1+rx nx]];
Print["At t = ",tf," approximate s = ",x0,",  error/dx^2 = ",(x0-x0f)/dx^2];
Print["At x = ",N[rx],", t = ",tf," approximate c = ",cx,",  error/dx^2 = ",(cx-u[rx,tf])/dx^2];

ListPlot[Table[{(i-1)dx,c[[i]]-u[(i-1)dx,tf]},{i,1,r}],AxesLabel->{"x","Error c[x,tf]"}]