%-----------------------------------------------------------------
% Script 'Discretization Methods', Example 3.9
% (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 = 1000;                         % 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 by forward difference method
% Temperature at time t=0 (initial condition)
T(1:n_s,1) = T_x0;

% Time integration, Eq. (3.135)
for k=1:1:n_t
  T(1,k+1) = T_0t;                  % bound. cond. for x=0
  T(n_s,k+1) = T_lt;                % bound. cond. for x=l
  for i=2:1:n_s-1                   % time int. step performed for each node
    T(i,k+1) = r*T(i-1,k)+(1.0-2.0*r)*T(i,k)+r*T(i+1,k);
  end
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')