D=qr1(A,epsilon)

function D=qr1(A,epsilon)

%Input - A is a symmetric tridiagona nxn matrix
%      - epsilon is the tolerance
%Output - D is the vector of eigenvalues

% NUMERICAL METHODS: MATLAB Programs
%(c) 1999 by John H. Mathews and Kurtis D. Fink
%To accompany the textbook:
%NUMERICAL METHODS Using MATLAB,
%by John H. Mathews and Kurtis D. Fink
%ISBN 0-13-270042-5, (c) 1999
%PRENTICE HALL, INC.
%Upper Saddle River, NJ 07458

%Initialize parameters

[n,n]=size(A);
m=n;

while (m>1)
   S=A(m-1:m,m-1:m);
   if abs(S(2,1))<epsilon
      A(m,m-1)=0;
      A(m-1,m)=0;
    
   else
      shift=eig(S);
      [j,k]=min([abs(A(m,m)-shift(1)) abs(A(m,m)-shift(2))]);
   end
   [Q,U]=qr(A-shift(k)*eye(n));
   A=U*Q+shift(k)*eye(n);
   m=m-1;
end

D=diag(A);