% This Matlab script calculates the eigenvalues of the one-dimensional
% convection-diffusion equation discretized by finite differences.
% The discretization uses central differences in space and forward
% Euler in time.
%
% Periodic bcs are set.
%

clear all;

% Set-up grid
xL = -4;
xR =  4;
Nx = 21; % number of points
x = linspace(xL,xR,Nx);

% Calculate cell size in control volumes (assumed equal)
dx = x(2) - x(1);

% Set velocity 
u = 1;

% Set viscosity using Reynolds based on dx
Re = 1e+2; % Don't set to zero identically... use a small positive value
mu = u*dx/Re;

% Set timestep
CFL = 0.5;
dt = CFL*dx/abs(u);

% Allocate matrix to hold stiffness matrix (A).
%
A = zeros(Nx-1,Nx-1);

% Construct A except for first and last row
for i = 2:Nx-2,
  A(i,i-1) =  u/(2*dx) + mu/dx^2;
  A(i,i  ) =          -2*mu/dx^2;
  A(i,i+1) = -u/(2*dx) + mu/dx^2;
end

% Periodic bcs
  
A(1   ,Nx-1) =  u/(2*dx) + mu/dx^2;
A(1   ,1   ) =          -2*mu/dx^2;
A(1   ,2   ) = -u/(2*dx) + mu/dx^2;
A(Nx-1,Nx-2) =  u/(2*dx) + mu/dx^2;
A(Nx-1,Nx-1) =          -2*mu/dx^2;
A(Nx-1,1   ) = -u/(2*dx) + mu/dx^2;
  

% Calculate eigenvalues of A
lambda = eig(A);

% Plot lambda*dt
plot(lambda*dt,'r*');
xlabel('Real \lambda\Delta t');
ylabel('Imag \lambda\Delta t');

% Overlay Forward Euler stability region
th = linspace(0,2*pi,101);
hold on;
plot(-1 + sin(th),cos(th));
hold off;
axis('equal');
grid on;