clear
Nk = 100;
k = linspace(0,pi,Nk);  % frequencies from DC to Nyquist (grid spacing = 1).

% Working in dimensionless units:
%     grid spacing Delta = 1
%     nominal phase speed c = 1

safety = 2.5         % stability safety factor
delta = 1/safety     % dimensionless timestep

% epsilon = delta^2 + 4*delta^4  % half-step bias
epsilon = (1 - 2*delta^2 - sqrt(1 - 4*delta^2)) / (2*delta^2) 
   % (exact stability threshold)

% Diffusion to mimic the RK2+ numerical diffusion.
D = (1+epsilon) * delta/10

% Eigenvalues:

E_pde = exp(i*k*delta);

E_ode = exp( i*delta*sqrt( 2*( 1 - cos(k) ) ) );

E_rk2 = 1 + i*delta*sqrt(2*(1 - cos(k))) - delta^2*(1 - cos(k));

E_rk2plus = 1 + i*delta*sqrt(2*(1 - cos(k))) - (1 + epsilon)*delta^2*(1 - cos(k));

E_diffusion = E_pde .* exp( -( k.^2 * D * delta / 2 ) );

% Graphics...
% myfont = '/sw/lib/X11/fonts/applettf/Arial.ttf'; 
   % Name of a font file that works for print(). 

figure(1)
plot(k, abs(E_pde),       'k-', 'linewidth', 2,...
     k, abs(E_diffusion), 'g-', 'linewidth', 2,...
     k, abs(E_ode),       'ob', 'linewidth', 2,...
     k, abs(E_rk2),       'r-', 'linewidth', 2,...
     k, abs(E_rk2plus),   'b-', 'linewidth', 2);
% axis([0,pi,0,1.1])
title('Eigenvalue Amplitude')
ylabel('Amplitude |{\xi}|')
xlabel('Wavenumber {k}')
legend('PDE','diffusion','Centered Differencing','RK2','RK2+ (stabilized)');
legend('location','northwest')
print('amplitude.pdf');

phi_pde = atan( imag(E_pde) ./ real(E_pde) );
phi_diffusion = atan( imag(E_diffusion) ./ real(E_diffusion) );
phi_ode = atan( imag(E_ode) ./ real(E_ode) );
phi_rk2 = atan( imag(E_rk2) ./ real(E_rk2) );
phi_rk2plus = atan( imag(E_rk2plus) ./ real(E_rk2plus) );

figure(2)
plot(k, phi_pde,       'k-', 'linewidth', 2,...
     k, phi_diffusion, 'g-', 'linewidth', 2,...
     k, phi_ode,       'ob', 'linewidth', 2,...
     k, phi_rk2,       'r-', 'linewidth', 2,...
     k, phi_rk2plus,   'b-', 'linewidth', 2);
title('Eigenvalue Phase')
ylabel('phase {\phi}')
xlabel('wavenumber {k}')
legend('PDE','diffusion','Centered Differencing','RK2','RK2+ (stabilized)');
legend('location','northwest')
print('phase.pdf');
