%-----------------------------------------------------------------
% Script 'Discretization Methods', Example 3.2
% (c) 09/2004 by A. Schmidt, IAM, Uni Stuttgart
%-----------------------------------------------------------------
clear all

% Variables
l = 1.0;                            % length of the rod
T = 1.0;                            % temperature of the basin
k = 1.0;                            % thermal conductivity of the rod
q = 1.0;                            % heat flux into the rod
n = 6;                              % number of nodes (without ghost points)

% Parameters
h = l/(n-1.0);                      % nodal distance

%-----------------------------------------------------------------
% Calculation of the System Matrix A: Az=b [Eq. (3.31)]
% first row (i=1)
A(1,1) = -2.0;
A(1,2) = 2.0;

% middle block
for i=2:1:n-2
  A(i,i-1) = 1.0;
  A(i,i) = -2.0;
  A(i,i+1) = 1.0;
end

% last row (i=n-1)
A(n-1,n-2) = 1.0;
A(n-1,n-1) = -2.0;

% Calculation of the right-hand side b: Az=b [Eq. (3.31)]
b(1) = -2.0*h/k*q;
b(n-1) = -T;
b = b';                             % make b a column vector

%-----------------------------------------------------------------
% Calculation of the vector of unknowns z (temperature): Az=b
z = A\b;

%-----------------------------------------------------------------
% Output
% sample points
for i=1:1:n
  x(i) = (i-1)*h;
end

% add temperature of the rod at the basin
z(n) = T;

% plot
figure(1)
plot(x,z,'k')
xlabel('Location / m')
ylabel('Temperature / °C')