%
%  lect22.m -- stability analysis
%
%- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 

clf

%  stab = c dt/dx

for stab = 0.5:0.5:1.5

%  basic parameters

N = 128;  L = pi;   	dx = 2*L/N;
c = 1.0;  			dt = stab*(dx/c);

%  wavenumber vector

n1 = -N/2 :1: N/2;

%  amplification factor

zu = 1 - (c*dt/dx)*(1 - exp(-i*(pi*dx/L)*n1));
zd = 1 + (c*dt/dx)*(1 - exp( i*(pi*dx/L)*n1));

%- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
%
%  plot here

subplot(2,1,1)
plot(n1/N,abs(zu),'b')
hold on

subplot(2,1,2)
plot(n1/N,abs(zd),'b')
hold on

end

subplot(2,1,1)
title('\bf stability of upwind differencing')
xlabel('\bf wavenumber (n/N)')
ylabel('\bf amplification (|Z|)')
axis([0 0.5 0 2])

subplot(2,1,2)
title('\bf instability of downwind differencing')
xlabel('\bf wavenumber (n/N)')
ylabel('\bf amplification (|Z|)')
axis([0 0.5 0 2])