The solution for a time dependent 2-D problem is computed for a set of grid points that cover the problem domain, a representative grid is shown in the diagram
The model for time dependent 2-D heat transfer is
Note that temperature in the plate is a function of three variables, T(t, x, y), and taking samples at discrete times and grid points gives values for T(n, i, j). The standard notation for the sampled values is Tn(i,j).
We make the usual Euler's substitution for the time derivative
The FDM substitutions for the spatial derivatives are
Making both substitutions into the differential equation model gives
and we can rearrage the equation to solve for Tn+1(i,j)
and we're good to go.
With dx = dy and C = dt*k / (cp*p*dx2), the code is
for n = 1:N
Tn = T; % save Tn before calculating Tn+1
for i = 2:sizeX - 1 % note that boundary points are omitted
for j = 2:sizeY - 1
T(i,j) = T(i,j) + C*(Tn(i+1,j) + Tn(i-1,j) + Tn(i,j+1) + Tn(i,j-1) - 4*T(i,j));
% apply boundary conditions here
Simple as that !