%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                 %
% Heat distribution on rectangular metal plate with constant      %
% temperature upper and lower boundaries (Dirichlet B.C.) and     %
% constant flows through the metal plate (Neumann B.C. on right   %
% and left) with Qx=32 and h=1.                                   %
%
%                       60                                        
%               
%       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)  
%          |                               |
% 32 -> T(3,1)  T(3,2)  T(3,3)  T(3,4)  T(3,5)  -> 32 
%          |                               |
%       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)  
%                                                                
%                       20                                       
%                                                                
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function gs_insulated

% create grid with boundary

T = zeros(5,5);

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

error

figure;        % plot temperature distribution
surf( T );
view(130,15)
title('Temperature distribution');

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