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

%-----------------------------------------------------------------
% remarks: - in contrast to the script, the temperature is plotted
%            against location AND time 
%          - keeping the given values, the critical time step
%            occurs for n_t = 200; check out n_t = 196 e.g.
%-----------------------------------------------------------------
clear all

% Variables
l = 1.0;                            % length of the rod
T_x0 = 0.0;                         % initial cond.: temperature T(x,0) 
T_0t = 100.0;                       % bound. cond.: temperature T(0,t) 
T_lt = 0.0;                         % bound. cond.: temperature T(l,t) 
alpha = 1.0e-5;                     % thermal diffusivity 
t_max = 1000.0;                     % time under consideration
n_s = 101;                          % number of spatial nodes
n_t = 5;                            % number of temporal nodes (time int.) 

% Parameters
dx = l/(n_s-1.0);                   % spatial nodal distance
dt = t_max/n_t;                     % time step size
r = alpha*dt/dx^2;                  % parameter r

%-----------------------------------------------------------------
% Calculation / time integration implicit by backward difference method
% Temperature at time t=0 (initial condition)
T(1:n_s,1) = T_x0;

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

% middle block
for i=2:1:n_s-3
  A(i,i-1) = -r;
  A(i,i) = 1.0+2.0*r;
  A(i,i+1) = -r;
end

% last row (i=n_s-1)
A(n_s-2,n_s-3) = -r;
A(n_s-2,n_s-2) = 1.0+2.0*r;

% Calculation of the right-hand side b2: Az=b=b1+b2
b2(1) = r*T_0t;
b2(2:n_s-3) = 0.0;
b2(n_s-2) = r*T_lt;

% Time integration, Eq. (3.157)
for k=1:1:n_t
  b1(1:n_s-2) = T(2:n_s-1,k);       % right-hand side b1
  b = (b1+b2)';                     % make b a column vector
  
  z = A\b;                          % Calculation of the vector of unknowns z 
                                    % (temperature): Az=b
  T(1,k+1) = T_0t;                  % bound. cond. for x=0
  T(n_s,k+1) = T_lt;                % bound. cond. for x=l
  T(2:n_s-1,k+1) = z;               % add solution to temperature array
end 

%-----------------------------------------------------------------
% Output
% sample points
for i=1:1:n_s
  x(i) = (i-1)*dx;
end
for i=1:1:n_t+1
  t(i) = (i-1)*dt;
end

% for a better plot transpose matrix T (swap t and x)
T=T';

% plot
figure(1)
mesh(x,t,T)
xlabel('x / m')
ylabel('time / s')
zlabel('Temperature / °C')