(* 
(1-0.1 Sin[π x])(((∂^2)u)/(∂(x^2)))=((∂u)/(∂t));
(1-0.1 u)(((∂^2)u)/(∂(x^2)))=((∂u)/(∂t));

u(x,0)=f[x];
u(0,t)=ϕ[t], u(1,t)=ψ(t) *)

(* a[x_]:=1.0-0.5Sin[π x] *)
(* a[u_]:=1- u *)
(* a[t_]:=1.0+ 1Sin[50 π t] *)
a[t_,u_]:=1.0+4 u Sin[50 π t]
ϕ[t_]:=0;
ψ[t_]:=0;
f[x_]:=-x(x-1);
tf=0.2;
n=20;
Δx=1.0/n;
eqs=Table[
D[Subscript[u,r][t],{t,1}]==a[t,Subscript[u,r][t]](Subscript[u,r+1][t]-2Subscript[u,r][t]+Subscript[u,r-1][t])/((Δx)^2),{r,1,n-1}]/.{Subscript[u,0][t]->ϕ[t],Subscript[u,n][t]->ψ[t]};
coi=Table[Subscript[u,r][0]==f[r Δx],{r,1,n-1}];
vars=Table[Subscript[u,r][t],{r,1,n-1}];
difeqs=Join[eqs,coi];

difeqs[[3]]

sol=NDSolve[difeqs,vars,{t,0,tf}];

nt=19;
Δt=tf/nt;
data={};
Do[
vx=r Δx;
vt=s Δt;
vu=Subscript[u,r][t]/.sol[[1]]/.{t->vt};
point={vx,vt,vu};
AppendTo[data,point];
,{r,1,n-1},{s,0,nt}]
Do[
vx=0;
vt=s Δt;
vu=ϕ[t];
point={vx,vt,vu};
AppendTo[data,point];
vx=1;
vt=s Δt;
vu=ψ[t];
point={vx,vt,vu};
AppendTo[data,point];
,{s,0,nt}];

ListPlot3D[data,PlotRange->All]

ListDensityPlot[data,PlotRange->All]

MatrixForm[data]