% waveeqn.m

clear all; close all;

% set various values
ell = 10 ;       % length of string
a=0.28  ;        % the left edge of the spike is at a*ell
d=0.3 ;          % the peak of the spike is at d*ell
b=0.32 ;         % the right edge of the spike is at b*ell
c=1.0 ;          % wave equation c (speed of propagation)
dx = 0.05 ;      % spacing in x space for plotting
dt = 0.05 ;      % time step for animation
tfinal = 20 ;    % final time for animation
nModes = 400 ;   % number of Fourier modes
delay = 0.0 ;    % extra time (in seconds) between animation frames 

% compute the Fourier coefficients for the initial position of the string

alpha = zeros(1, nModes) ;
for k = 1:nModes
     alpha(k) = - 2*sin(k*pi*a)/( (d-a)*k*k*pi*pi ) ...
                + 2*(b-a)*sin(k*pi*d)/( (d-a)*(b-d)*k*k*pi*pi ) ...
                - 2*sin(k*pi*b)/( (b-d)*k*k*pi*pi ) ;
end

% compute the temporal frequencies

for k = 1:nModes
     Om(k) = c*k*pi/ell ;
end

% run the time evolution

X = 0:dx:ell ;
set(gcf,'DoubleBuffer', 'on') ;
for t = 0:dt:tfinal                     % for each fixed t
     B = alpha .* cos(t*Om) ;           % construct coeffs of sin(k*x*pi/ell)
                                        % for time t
     U = zeros(1, length(X)) ;          % next construct u(x,t)
     for k = 1:nModes 
          U = U + B(k)*sin( k*pi*X/ell ) ;
     end
     plot(X,U) ;                        % plot u(x,t)
     axis([0 ell -1 1]) ;
     title(sprintf('t = %4.1f', t)) ; 
     pause(delay) ;
end