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

% Variables
l = 0.05;                           % length of the rod
R = 0.01;                           % radius of the rod
T_wall = 330.0;                     % temperature of the wall
T_oo = 30.0;                        % ambient temperature
k = 50.0;                           % thermal conductivity of the rod
h = 100.0;                          % convection coefficient (rod-air)
n = 6;                              % number of nodes (without ghost points)
n_a = 200;                          % number of sample points 
                                    % (analytical solution)

% Parameters
dx = l/(n-1.0);                     % nodal distance
m = sqrt(2.0*h/(R*k));              % factor m

%-----------------------------------------------------------------
% Calculation of the System Matrix A: Az=b [Eqs. (3.58), (3.61)]
% first row (i=2)
A(1,1) = -(2.0+m^2*dx^2);
A(1,2) = 1.0;

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

% last row (i=n)
A(n-1,n-2) = 2.0;
A(n-1,n-1) = -(2.0+m^2*dx^2);

% Calculation of the right-hand side b: Az=b [Eq. (3.54)]
b(1) = -(T_wall-T_oo);
b(2:n-1) = 0.0;
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) = (T_wall-T_oo)*cosh(m*(l-x_a(i)))/cosh(m*l);
end

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

% adjusted vector of (absolute) temperatures
z = z+T_oo;                         % add T_oo (absolute value)
z_a = z_a+T_oo;                     % also for analytical solution

for i=n:-1:2                        % sample point at wall has to be added
  z(i) = z(i-1);
end
z(1) = T_wall;

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