function [Lambda,A1,A2]=qr96v(A,Nmax,TOL);
%%function [Lambda,A1,A2]=qr96v(A,Nmax,TOL);
%%   The QR algorithm is applied to a symmetric tridiagonal matrix A to
%%   obtain either the eigenvalues of A, or some splitting of A
%%   This is verbose mode, the decomposition is shown at each stage.
%%  INPUT
%%    A - a symmetric tridiagonal matrix
%%    Nmax - maximum number of iterations
%%    TOL - tolerance
%%  OUTPUT
%%    Lambda - column of approximate eigenvalues
%%      A1  -  kxk symmetric tridiagonal matrix
%%      A2  -  (n-k)x(n-k) symmetric tridiagonal matrix
%%    The algorithm may be applied once again to A1 and A2 separately

%% Step 1
k=1;shift=0;[n,m]=size(A);FLAG=0;
%% Step 2
while (k<=Nmax & FLAG==0)
    for j=2:n;
       b(j)=abs(A(j,j-1));
    end;
    if b(n)<=TOL;
        Lambda=[Lambda;A(n,n)+shift];
        A=A(1:n-1,1:n-1);n=n-1;
    end;
    I=find(b(3:n)<=TOL);
    %if isempty(I)==0; %%Use with Matlab <=3.5
     if size(I)~=[0,0];
        FLAG=1;
        A1=A(1:I(1)+1,1:I(1)+1);A2=A(I(1)+2:n,I(1)+2:n);    
    end;
    if b(2)<= TOL;
          Lambda=[Lambda;A(1,1)+shift];
          A=A(2:n,2:n); n=n-1;
    end;
    %Step 4
    b=-(A(n-1,n-1)+A(n,n));c=A(n,n)*A(n-1,n-1)-(A(n-1,n)).^2;
    d=sqrt(b^2-4*c);
    if b>0;
       mu1=-2*c/(b+d);mu2=-(b+d)/2;
    else
       mu1=(d-b)/2;mu2=2*c/(d-b);
    end;
    if n==2;
      Lambda=[Lambda;mu1+shift;mu2+shift];
      FLAG=1;
    end;
    s=min([abs(mu1-A(n,n));abs(mu2-A(n,n))])+A(n,n);
    shift=shift+s; %% Step 5
fprintf(' stage  is \n ')
k
fprintf('shift  is \n ')
s
    %% Step 6   
    AA=A-s*eye(n); 
    %% Step 7
%%% Here we do extra work  to show R and Q
    x(1)=AA(1,1);y(1)=AA(2,1);
    for j=2:n;
        z(j-1)= sqrt(x(j-1).^2+(AA(j,j-1)).^2);
        c(j)=x(j-1)/z(j-1);
        s(j)=AA(j,j-1)/z(j-1);
        p=j-1;
        if p==1
            Q=[[c(2) s(2);-s(2) c(2)] zeros(2,n-2);zeros(n-2,2) eye(n-2)]';
        elseif p<n-1
            Q=Q*[eye(p-1) zeros(p-1,2) zeros(p-1,n-p-1);
              zeros(2,p-1) [c(p+1) s(p+1);-s(p+1) c(p+1)] zeros(2,n-p-1);
              zeros(n-p-1,p-1) zeros(n-p-1,2) eye(n-p-1)]';
        elseif p==n-1
            Q=Q*[eye(p-1) zeros(p-1,2) ;
              zeros(2,p-1) [c(p+1) s(p+1);-s(p+1) c(p+1)]]';
        end;
        q(j-1)=c(j)*y(j-1)+s(j)*AA(j,j);
        x(j)= -s(j)*y(j-1)+c(j)*AA(j,j);
        if j<n
            r(j-1)=s(j)*AA(j+1,j);
            y(j)=c(j)*AA(j+1,j);
        end;
    end;
fprintf('Q for A-sI  at this stage is\n')
Q
pause
for j=1:n-2;
    R(j,j)=z(j);R(j,j+1)=q(j);R(j,j+2)=r(j);
end;R(n-1,n-1)=z(n-1);R(n-1,n)=q(n-1);R(n,n)=x(n);
fprintf('R for A-sI  at this stage is\n')
R
pause
clear Q U R
    %% Step 8
    z(n)=x(n);
    A(1,1)=s(2).*q(1)+c(2).*z(1);
    A(2,1)=s(2)*z(2);A(1,2)=A(2,1);
    for j=2:n-1;
      A(j,j)=s(j+1).*q(j)+c(j).*c(j+1).*z(j);
      A(j+1,j)=s(j+1).*z(j+1);A(j,j+1)=A(j+1,j);
    end;
    A(n,n)=c(n)*z(n);
    k=k+1;
fprintf('New A at this stage is\n')
A
pause
end;
if (k==Nmax & FLAG==0)
  message = 'Maximum number of iteration exceeded without success'
end;
