% Set non-dimensional thermal coefficient
k   = 1.0; % this is really k/(rho*cp)

% Set length of bar
L = 1.0; % non-dimensional

% Set up grid size
Nx = 1000;
h = L/Nx;
x = linspace(0,L,Nx+1);

% Initialize a sparse matrix to hold stiffness matrix
A = spalloc(Nx-1,Nx-1,3*(Nx-1));

% Calculate stiffness matrix

for ii = 1:Nx-1,
    
  if (ii > 1),
    A(ii,ii-1) = k/h^2;
  end
  
  if (ii < Nx-1),
    A(ii,ii+1) = k/h^2;
  end
  
  A(ii,ii) = -2*k/h^2;
  
end

% Find eigenvalues of A
lam = eig(A);

% Determine minimum and maximum eigenvalues of A
[lammin, imin] = min(abs(lam)); % Find minimum magnitude eigenvalue
[lammax, imax] = max(abs(lam)); % Find maximum magnitude eigenvalue

% Determine spectral condition number
scn = lammax/lammin;

% Print out info to screen
fprintf('Minimum eigenvalue = %f\n',lam(imin));
fprintf('Maximum eigenvalue = %f\n',lam(imax));
fprintf('Spectral condition number = %f\n',scn);

% Plot eigenvalues
plot(real(lam),imag(lam),'*');
xlabel('Real \lambda');
ylabel('Imag \lambda');