(* 
(((∂^2)u)/(∂(x^2)))+(((∂^2)u)/(∂(y^2)))=((∂u)/(∂t));

u(x,y,0)=f[x,y];
u(0,y,t)=ϕ[y,t], u(1,y,t)=ψ(y,t);
u(x,0,t)=α[x,t], u(x,1,t)=β(x,t);

(Subscript[u, rx,ry])(t)=u(rx Δx,ry Δy,t)


 *)

Plot3D[Sin[π x]Sin[π y],{x,0,1},{y,0,1}]

ϕ[y_,t_]:=0;
ψ[y_,t_]:=0;
f[x_,y_]:=Sin[π x]Sin[π y];
tf=0.2;
nx=4;
ny=4;
Δx=1.0/nx;
Δy=1.0/ny;
eqs=Flatten[Table[
D[Subscript[u,rx,ry][t],{t,1}]==(Subscript[u,rx+1,ry][t]-2Subscript[u,rx,ry][t]+Subscript[u,rx-1,ry][t])/((Δx)^2)+(Subscript[u,rx,ry+1][t]-2Subscript[u,rx,ry][t]+Subscript[u,rx,ry-1][t])/((Δy)^2),{rx,1,nx-1},{ry,1,ny-1}]];

coi=Flatten[Table[Subscript[u,rx,ry][0]==f[rx Δx,ry Δy],{rx,1,nx-1},{ry,1,ny-1}]];

vars=Flatten[Table[Subscript[u,rx,ry][t],{rx,1,nx-1},{ry,1,ny-1}]]

difeqs=Join[eqs,coi]

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]