function [D,hist]=qr2_enh(A,epsilon,maxit)

% PROGRAM FROM: MATHEWS FINK, P. 606.

% INPUT: A IS A SYMMETRIC TRIDIAGONAL n x n MATRIX.
%        epsilon IS THE TOLERANCE.
% OUTPUT: D IS THE n x 1 VECTOR OF EIGENVALUES.
%         hist=ITERATIONS VECTOR:

% INITIALIZE PARAMETERS.
[n,n]=size(A);
m=n;
D=zeros(n,1);
B=A;

index=1;

while ( m > 1 )
    
    hist(index,1)=0; hist(index,2)=0;
    
    while (abs(B(m,m-1)) >= epsilon) & (hist(index) < maxit)
        
        % CALCULATE SHIFT.
        S=eig(B(m-1:m,m-1:m));
        [j,k]=min([abs(B(m,m))]*[1 1]'-S);
        
        % QR FACTORIZATION OF B.
        [Q,U]=qr(B-S(k)*eye(m));
        
        % CALCULATE NEXT B.
        B=U*Q+S(k)*eye(m);
        
        hist(index,1)=hist(index,1)+1;
        hist(index,2)=abs(B(m,m-1));
    end
    
    if (hist(index) == maxit)
        B
    end
    
    %PLACE mth EIGENVALUE IN A(m,n).
    A(1:m,1:m)=B;
    
    % REPEAT PROCESS ON THE m-1 x m-1 SUBMATRIX OF A.
    m=m-1;
    B=A(1:m,1:m);
    
    index=index+1;
    
end

D=diag(A);