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

%-----------------------------------------------------------------
% remarks: - forces f=0 along the rod except for the free end
%          - in contrast to the figure in the script, the 
%            direction of the x-axis is opposed
%          - cross section area is denoted by 'a' since 'A' is the 
%            system matrix (Az=b)
%-----------------------------------------------------------------
clear all

% Variables
l = 1.0;                            % length of the rod
a = 1.0e-4;                         % cross section area of the rod
E = 2.1e+11;                        % Young's modulus
rho = 7.85e+3;                      % density of the material
N = 1.0e+3;                         % external load
dt = 1.0e-5;                        % time step size
t_max = 1.0e-3;                     % time under consideration
n_s = 101;                          % number of spatial nodes
u_x0 = 0.0;                         % initial cond.: displacement u(x,0) 
v_x0 = 0.0;                         % initial cond.: velocity v(x,0) 
u_0t = 0.0;                         % bound. cond.: displacement u(0,t) 

% Parameters
n_t = floor(t_max/dt);              % number of temporal nodes (time int.) 
dx = l/(n_s-1.0);                   % spatial nodal distance
r2 = E/rho*dt^2/dx^2;               % parameter r^2

%-----------------------------------------------------------------
% Calculation / time integration using an implicit scheme
% Displacement at time t=0 (initial condition)
u(1:n_s,1) = u_x0;

% ghost points displacements of all nodes Eq. (3.231)
u_gp = u_x0-dt*v_x0;

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

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

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

b(1:n_s-1) = 0.0;                   % initalize right-hand side b ...
b = b';                             % ... and make b a column vector

% Time integration
% first time step using ghost points (k=1)
u(1,2) = u_0t;                      % bound. cond. (displacement) for x=0

b(1) = 2.0*u(2,1)-u_gp+r2*u_0t;
b(2:n_s-2) = 2.0*u(3:n_s-1,1)-u_gp;
b(n_s-1) = 2.0*u(n_s,1)-u_gp+r2*2.0*N*dx/E/a;

z = A\b;                            % Calculation of the vector of unknowns z
u(2:n_s,2) = z;                     % add solution to displacement array

% time stepping (k>1)
for k=2:1:n_t
  u(1,k+1) = u_0t;                  % bound. cond. (displacement) for x=0
  
  b(1) = 2.0*u(2,k)-u(2,k-1)+r2*u(1,k+1);        % right-hand side b
  b(2:n_s-2) = 2.0*u(3:n_s-1,k)-u(3:n_s-1,k-1);
  b(n_s-1) = 2.0*u(n_s,k)-u(n_s,k-1)+r2*2.0*N*dx/E/a;

  z = A\b;                          % Calculation of the vector of unknowns z
  u(2:n_s,k+1) = z;                 % add solution to displacement 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

% plot
figure(1)
mesh(t,x,u)
xlabel('time / s')
ylabel('x / m')
zlabel('displacement / m')