(*mesh parameters*)
domain={{0,5},{0,5}};
α=.1;
Δt=.05; (*time step parameter*)
Δx=.5; (*both x and y direction mesh parameter*)


λ=((α*Δt)/((Δx)^2));
(*generate node records*)
xcoords=Table[domain[[1,1]]+i*Δx,{i,0,(domain[[1,2]]-domain[[1,1]])/Δx}];
ycoords=Table[domain[[2,1]]+i*Δx,{i,0,(domain[[2,2]]-domain[[2,1]])/Δx}];
nodeRecord=Table[{i,
xcoords[[Floor[(i-1)/Length[ycoords]]+1]],
ycoords[[Mod[i,Length[ycoords],1]]]},
{i,1,Length[xcoords]*Length[ycoords]}];

(*load the matrix*)
A=(1-4λ)*IdentityMatrix[Length[xcoords]*Length[ycoords]];
(*corners first*)
(*bottom left*)
A[[1,2]]=λ;
A[[1,1+Length[ycoords]]]=λ;
(*top left*)
A[[Length[ycoords],Length[ycoords]-1]]=λ;
A[[Length[ycoords],2*Length[ycoords]]]=λ;
(*bottom right*)
A[[(Length[xcoords]-1)*Length[ycoords]+1,(Length[xcoords]-2)*Length[ycoords]+1]]=λ;
A[[(Length[xcoords]-1)*Length[ycoords]+1,(Length[xcoords]-1)*Length[ycoords]+2]]=λ;
(*top right*)
A[[(Length[xcoords]-1)*Length[ycoords]+Length[ycoords],(Length[xcoords]-1)*Length[ycoords]]]=λ;
A[[(Length[xcoords]-1)*Length[ycoords]+Length[ycoords],(Length[xcoords]-1)*Length[ycoords]+Length[ycoords]-1]]=λ;
(*sides*)
(*left side*)
For[i=2,i<Length[ycoords],i++,
A[[i,i-1]]=λ;
A[[i,i+1]]=λ;
A[[i,i+Length[ycoords]]]=λ;
]
(*bottom side*)
For[i=2,i<Length[xcoords],i++,
A[[(i-1)*Length[ycoords]+1,(i-1)*Length[ycoords]+2]]=λ;
A[[(i-1)*Length[ycoords]+1,(i-2)*Length[ycoords]+1]]=λ;
A[[(i-1)*Length[ycoords]+1,i*Length[ycoords]+1]]=λ;
]
(*top side*)
For[i=2,i<Length[xcoords],i++,
A[[(i-1)*Length[ycoords]+Length[ycoords],(i-1)*Length[ycoords]+Length[ycoords]-1]]=λ;
A[[(i-1)*Length[ycoords]+Length[ycoords],(i-2)*Length[ycoords]+Length[ycoords]]]=λ;
A[[(i-1)*Length[ycoords]+Length[ycoords],i*Length[ycoords]+Length[ycoords]]]=λ;
]
(*right side*)
For[i=2,i<Length[ycoords],i++,
A[[i+(Length[xcoords]-1)*Length[ycoords],i+(Length[xcoords]-1)*Length[ycoords]-1]]=λ;
A[[i+(Length[xcoords]-1)*Length[ycoords],i+(Length[xcoords]-1)*Length[ycoords]+1]]=λ;
A[[i+(Length[xcoords]-1)*Length[ycoords],i+(Length[xcoords]-1)*Length[ycoords]-Length[ycoords]]]=λ;
]
(*middle*)
For[i=2,i<Length[ycoords],i++,
For[j=2,j<Length[xcoords],j++,
A[[Length[ycoords]*(j-1)+i,Length[ycoords]*(j-1)+i+1]]=λ;
A[[Length[ycoords]*(j-1)+i,Length[ycoords]*(j-1)+i-1]]=λ;
A[[Length[ycoords]*(j-1)+i,Length[ycoords]*j+i]]=λ;
A[[Length[ycoords]*(j-1)+i,Length[ycoords]*(j-2)+i]]=λ;
]
]

(*boundary conditions go here*)
(*leave blank for full conductivity outwards*)
(*ex. insulate all sides*)
For[i=1,i≤Length[A],i++,
A[[i,i]]+=1-Sum[A[[i,k]],{k,1,Length[A]}];
]

(*ex. left hand side is open to an infinite heat sink*)
For[i=1,i≤Length[ycoords],i++,
A[[i]]=ConstantArray[0,Length[A]];
]
(*end boundary conditions*)

(*initial conditions*)
values=Table[Append[Delete[nodeRecord[[i]],1],RandomReal[1]],{i,1,Length[xcoords]*Length[ycoords]}];
ListPlot3D[values,ColorFunction->"TemperatureMap",ColorFunctionScaling->False,Mesh->None]

(*visualization 1:*)
(*generate table of plots for successive iterations*)
showIterations={1,30};(*show iterations from {start,finish}*)
showEvery=5;(*show every 5 iterations*)
record={values};
iteratedvalues=values;
For[i=1,i≤showIterations[[2]],i++,
temp=Transpose[iteratedvalues];
u=temp[[3]];
u=A.u;
temp[[3]]=u;
iteratedvalues=Transpose[temp];
record=Append[record,iteratedvalues];
]
Table[ListPlot3D[record[[i]],ColorFunction->"TemperatureMap",ColorFunctionScaling->False,Mesh->None],{i,showIterations[[1]],showIterations[[2]],showEvery}]

(*visualization 2:*)
(*animation: run cell, click on the plot to stop animation*)
record={values};
iteratedvalues=values;
DynamicModule[{go=True},
Monitor[While[go,
temp=Transpose[iteratedvalues];
u=temp[[3]];
u=A.u;
temp[[3]]=u;
iteratedvalues=Transpose[temp];
record=Append[record,iteratedvalues];
p=ListPlot3D[iteratedvalues,ColorFunction->"TemperatureMap",ColorFunctionScaling->False,Mesh->None,PlotRange->Append[domain,{0,1}]];
],EventHandler[
p,{"MouseClicked":>(go=False)}]]
];
p