  function ShowTriDiagEigs()
% Illustrate the 3-term recurrence that can be used to evaluate the
% characteristic polynomial of a symmetric tridiagonal matrix.
close all
% Generate a random symmetric tridiagonal matrix...
n = 5;
A = randn(n,n); A = A+A'; A = tril(triu(A,-1),1);
% Extract the diagonal and subdiagonal...
a = diag(A);
b = diag(A,-1);
% Plot the characteristic polynomial of A
N = 200;
x = linspace(-norm(A)-.1,norm(A)+.1,N);
pVals = zeros(1,N);
for j=1:N
    pVals(j) = Eval_p(x(j),a,b);
end
plot(x,pVals,[x(1) x(N)],[0 0],'k',eig(A),zeros(1,n),'om')
title('Characteristic polynomial of a symmetric tridiagonal matrix') 
shg

  function z = Eval_p(lambda,a,b)
% Let A be a symmetric tridiagonal matrix with diagonal a(1:n) and
% subdiagonal b(1:n-1). Let lambda be a scalar. z = det(A - lambda*I).

% Idea. If pk(x) = det(A(1:k,1:k) - x*eye(k,k)), then
%             p_{0}    = 1         (convention)
%             p_{1}(x) = a(1) - x
%             p_{k}(x) = (a(k)-x)*p_{k-1}(x) - b(k-1)^2 p_{k-2}(x)  k=2:n
n = length(a);
pkm1 = 1;
pk = a(1) - lambda;
for k=2:n
    pkm2 = pkm1;
    pkm1 = pk;
    pk = (a(k) - lambda)*pkm1 - b(k-1)^2*pkm2;
end
z = pk;
    