(* 
(((∂^2)u)/(∂(x^2)))=((∂u)/(∂t));
u(x,0)=f[x];
u(0,t)=ϕ[t], u(1,t)=ψ(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}]==(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}];

eqs

coi

difeqs=Join[eqs,coi]

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

Plot[Subscript[u,3][t]/.sol[[1]],{t,0,tf},PlotRange->All]

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

MatrixForm[data]

ListPlot3D[data,PlotRange->All]

ListDensityPlot[data,PlotRange->All]