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

%-----------------------------------------------------------------
% remark: in contrast to the script, the displacement field is
%         only plottet for one point of time (t_max). All other
%         points of time are buffered in matrix u.
%-----------------------------------------------------------------
clear all

% Variables
l = 100.0;                          % length
c = 1.0;                            % wave speed (c>0)
dt = 0.5;                           % time step size
u_x0 = 0.0;                         % initial cond.: displacement u(x,0) 
u_0t = 1.0;                         % bound. cond.: displacement u(0,t) 
u_lt = 0.0;                         % bound. cond.: displacement u(l,t) 
t_max = 70.0;                       % time under consideration
n_s = 1001;                         % number of spatial nodes

% Parameters
n_t = floor(t_max/dt);              % number of temporal nodes (time int.) 
dx = l/(n_s-1.0);                   % spatial nodal distance
Co = c*dt/dx;                       % Courant number

%-----------------------------------------------------------------
% Calculation / time integration using Crank-Nicolson scheme
% Displacement at time t=0 (initial condition)
u(1,1) = u_0t;
u(2:n_s-1,1) = u_x0;
u(n_s,1) = u_lt;

% Calculation of the System Matrix A: Az=b [from Eq. (3.201)]
% first row (i=2)
A(1,1) = 4.0;
A(1,2) = Co;

% middle block
for i=2:1:n_s-3
  A(i,i-1) = -Co;
  A(i,i) = 4.0;
  A(i,i+1) = Co;
end

% last row (i=n_s-1)
A(n_s-2,n_s-3) = -Co;
A(n_s-2,n_s-2) = 4.0;

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

% Time integration
for k=1:1:n_t
  b(1) = 4.0*u(2,k)-Co*u(3,k)+2.0*Co*u_0t;        % right-hand side b
  b(2:n_s-3) = 4.0*u(3:n_s-2,k)+Co*(u(2:n_s-3,k)-u(4:n_s-1,k));
  b(n_s-2) = 4.0*u(n_s-1,k)+Co*u(n_s-2,k)-2.0*Co*u_lt;

  z = A\b;                          % Calculation of the vector of unknowns z 

  u(1,k+1) = u_0t;                  % bound. cond. (displacement) for x=0
  u(n_s,k+1) = u_lt;                % bound. cond. (displacement) for x=l
  u(2:n_s-1,k+1) = z;               % add solution to displacement array
end

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

% plot
figure(1)
plot(x,u(:,n_t+1))
xlabel('x / m')
ylabel('displacement / m')