%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Solution of Laplace (u_xx + u_yy = 0 ) on the following system:
% 
% Heat distribution on rectangular metal plate with
% constant temperature boundaries (Dirichlet)
%              
%                       100
%               
%      T(1,1)--T(1,2)--T(1,3)--T(1,4)--T(1,5)
%         |                               |
%      T(2,1)  T(2,2)  T(2,3)  T(2,4)  T(2,5)  
%         |                               |
%  75  T(3,1)  T(3,2)  T(3,3)  T(3,4)  T(3,5)  50  
%         |                               |
%      T(4,1)  T(4,2)  T(4,3)  T(4,4)  T(4,5) 
%         |                              |
%      T(5,1)--T(5,2)--T(5,3)--T(5,4)--T(5,5)  
%
%                        10

%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function gauss_seidel

% create grid with boundary

T = zeros(5,5);

T(:,1) = 75*ones(1,5);
T(:,5) = 50*ones(1,5);
T(5,:) = 10*ones(1,5);
T(1,:) = 100*ones(1,5);

T_new = zeros(5,5); % create new iteration matrix

% main iteration

error = 1;

while ( error > 0.001 )
  for i=2:4
    for j=2:4
      T_new(i,j) = ( T(i+1,j) + T(i-1,j) + T(i,j+1) + T(i,j-1) )/4;
    end
  end
  error = max(max( abs(T(2:4,2:4) - T_new(2:4,2:4)) ));
  T(2:4,2:4) = T_new(2:4,2:4);    % update T

  error
  
end

T

% plot temperature distribution

figure;       
surf( T );
view(-220,20)
title('Temperature distribution');

% plot isothermal lines

figure;       
contour(T')
title('Isothermal lines')