%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% 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 laplace

% create operator matrix m  
  
m = diag( -4*ones(9,1) );
band1 = [ 1 1 0 1 1 0 1 1 ];
m = m + diag( band1, -1 ) + diag( band1, 1 );
m = m + diag( ones(6,1), -3 ) + diag( ones(6,1), 3 ); % Dirichlet

% construct vector b using boundary conditions

bndry = [ -175 -100 -150 -75 0 -50 -85 -10 -60 ]';

% find solution using matrix inversion

res = inv(m) * bndry;

% convert solution from vector to matrix format

mres = [ res(1:3)'; res(4:6)'; res(7:9)' ];

% attach boundary

t_dist(:,1) = 75*ones(1,5);
t_dist(:,5) = 50*ones(1,5);
t_dist(1,:) = 100*ones(1,5);
t_dist(5,:) = 10*ones(1,5);
t_dist(2:4,2:4) = mres

% plot temperature distribution

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

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