% 2-D time dependent heat transfer
% model dT/dt = k*(d2T/dx2 + d2T/dy2)
% with Dirichlet and Neumann boundary conditions 
% and several heat sources/sinks

k = 2;

L = 20;

T = zeros(L, L);    % the estimated temperature values
To = zeros(L, L);    % old estimated temperature values 
f = zeros(L, L);    % heat source/sink values

dx = 1; dy = 1;

xi=0:dx:L; yi=0:dx:L;

dt = .1  % the time subinterval lenght
N = 2000;  % number of iterations

f(5,5) = 100;
f(15,10) = 100;

% initial temperature values at all grid points are 0 (may be overriden by Dirichlet BC)

% initialize Dirichlet boundary values
for i=[1 L]
    for j=1:L
        if i==1 T(i,j) = 0;
        else T(i,j) = 250;
        end
    end
end

contour(T)
axis([0 20 0 20])
caxis([0 300])
colorbar
keyboard


for n=1:N        % step through time
   % compute temp for all interior points
    To = T;
    for i=2:L-1
        for j=2:L-1
            T(i,j) = To(i,j) + (dt/dx^2)*k*(To(i+1,j) + To(i-1,j) + To(i,j+1) + To(i,j-1) - 4*To(i,j)) + f(i,j);
        end
    end

   % compute Neumann boundary values
    for j=[1 L]
        for i=2:L-1
            if j==1 T(i,j) = T(i,2);
            else T(i,j) = T(i,L-1);
            end
        end
    end
    
    if mod(n,20)==0
        %surf(T)
        %axis([0 20 0 20 0 300])
        contour(T)
        caxis([0 300])
        colorbar
        keyboard
    end
end