%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Solution of Laplace (u_xx + u_yy = 0 ) on the following system:
% 
% Heat distribution on rectangular metal plate with three
% constant temperature boundaries (Dirichlet) and one Neumann.
%              

%              
%                       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)  
%
%                     isolated
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function gs_insulated

% create grid with boundary
  
T = zeros(5,5);
T(1,:) = 100 * ones(1,5);
T(:,1) = 75 * ones(1,5);
T(:,5) = 50 * ones(1,5);

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

% main iteration

for n=1:60
  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
    T_new(5,i) = ( T(5,i+1) + T(5,i-1) + 2 * T(4,i) )/4;
  end
  error = abs(max(max( T(2:5,2:4) - T_new(2:5,2:4)) ));
  T(2:5,2:4) = T_new(2:5,2:4);    % update T
end

error

figure;        % plot temperature distribution
surf( T );
view(140,25)
title('Temperature distribution');

figure;        % plot isothermal lines
contour(T')
title('Isothermal lines')