%-----------------------------------------------------------------
% Script 'Discretization Methods', Example 3.12
% (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.05;                          % 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 = 80.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 by forward difference (time) / backward 
% difference method (space)
% Displacement at time t=0 (initial condition)
u(1:n_s,1) = u_x0;

% Time integration, Eq. (3.181)
for k=1:1:n_t
  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
  for i=2:1:n_s-1                   % time int. step performed for each node
    u(i,k+1) = u(i,k) - Co*(u(i,k)-u(i-1,k));
  end
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')