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
First analytic solution Ac is found such that 
s(t) = (Sqrt[1-130(t^3)]),  (* Example 3 *) 
Ac=c(x,t) = ((s(t))^2) (∑_(n=2)^∞)(((x/s(t)-1))^n)(Subscript[P, n])(t)/n!
coefficients  of (Subscript[P, n])(t)  stored in vectors (Subscript[p, n])
Then the numerical method is used to obtain c using 
given Ac(0,t), Ac(x,0)

nmax=100;
p=Table[0,nmax];
p[[2]]={1};
p[[3]]={0,0,195};
For[n=4,n≤nmax,n++,
nm=2n-3;pl=Table[0,nm];
For[i=3,i≤nm,i++,pl[[i]]=195p[[n-1]][[i-2]]];
For[i=3,i≤nm-2,i++,pl[[i]]+=195(n-4)p[[n-2]][[i-2]]];
For[i=2,i≤nm-4,i++,
pl[[i-1]]+=(i-1)p[[n-2]][[i]];
pl[[i+2]]-=130(i-1)p[[n-2]][[i]]];
p[[n]]=pl];

Print["p_5 = ",p[[5]]];

Ac[x_,t_]:=Module[{x02=1-130 t^3,x0,z,y,s,,b,er=1, eps=10.^(-13),nm,go=5,n=3,i},

x0=x02^.5;y=x/x0-1;
z=y y .5; s=z;
While[go>0,
z=y z/n;nm=2n-3;
b=p[[n]][[nm]];
For[i=nm-1,i>0,i--,b=b t+p[[n]][[i]]];
er=z  b;
s=s+er;n++;
If[Abs[er]<eps,go--,go=5; (* if |er|< eps in five consecutive times quit *)
If[n>nmax,Print["No convergence"];go=-1]]];
s x02];


Print["Analytic value of c(x,t) at x=.2, t=.16: ",Ac[.2,.16]];

Using the numerical method to solve 
((∂c)/(∂t))=(((∂^2)c)/(∂(x^2)))-1      for   0 < x < s(t),    0 < t < tf
((∂c)/(∂x))(s(t),t)=c(s(t),t)=0,   for 0 < t < tf
c(0,t)  = Ac(0,t)    for   0 < t < tf                    (Ac given above)
c(x,0) = Ac(x,0)   for  0 < x < 1

tf=.16;                (* final time *)

(* at tf to be compared with:*)
x0f=(1-130 tf^3)^.5;            (* exact value of s(tf) *)
rx=2/10;                          (* reporting at x=rx at tf, should be faction *)
vf=Ac[rx,tf];              (* exact value of c(rx,tf) *)


wa=0.5;
wb=1.5;
Print["Window [ ",wa," , ",wb," ]dx"];

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

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

c=Table[Ac[i dx,0],{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;
t+=dt;
cn[[1]]=Ac[0,t];    (* for prescribed boundary value at x=0 *)

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-130(i dt)^3)^.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-Ac[rx,tf])/dx^2];

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