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

% Variables
l = 1.0;                            % length of the rod
T0 = 0.0;                           % bound cond: temperature at x=0
Tl = 0.0;                           % bound cond: temperature at x=l
k = 1.0;                            % thermal conductivity of the rod
n = 6;                              % number of nodes (without ghost points)
n_a = 200;                          % number of sample points 
                                    % (analytical solution)

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

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

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

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

% Calculation of the right-hand side b: Az=b [Eq. (3.38)]
for i=1:1:n-2
  x(i) = i*h;
  b(i) = -h^2/k*sin(pi*x(i)/l);
end
b = b';                             % make b a column vector

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

%-----------------------------------------------------------------
% Analytical solution
for i=1:1:n_a
  x_a(i) = (i-1)*l/(n_a-1);         % sample point location
  z_a(i) = l^2/pi^2/k*sin(pi*x_a(i)/l);
end

%-----------------------------------------------------------------
% Output
for i=n-1:-1:2                      % sample points have to be added
  x(i) = x(i-1);
  z(i) = z(i-1);
end
x(1) = 0.0;
x(n) = l;
z(1) = T0;
z(n) = Tl;

% plot
figure(1)
plot(x,z,'ro',x_a,z_a,'k')
xlabel('Location / m')
ylabel('Temperature / °C')
legend('numerical','analytical')
title('o  numerical      --  analytical')