  function e = EigsInInterval(T,L,R,nL,nR)
% T is an unreduced symmetric tridiagonal matrix
% L < R and neither L nor R are eigenvalues of T
% nL is the number of eigenvalues to the left of L
% nR is the number of eigenvalues to the left of R
% e is a column vector (possibly empty) consisting of all the eigenvalues in 
% the interval [L,R] computed to full machine precision (as in the
% function Sturm).

if nargin ==3
    % Process the top level call...
    nL = EigsToLeft(T,L);
    nR = EigsToLeft(T,R);
end
if nL==nR
    % No eigenvalues in the interval..
    e = [];
elseif nL == nR-1
    % One eigenvalue in the interval...
    e = Sturm(T,nR,L,R);
else
    % Halve the interval and recur..
    Mid = (L+R)/2;
    nMid = EigsToLeft(T,Mid);
    eL = EigsInInterval(T,L,Mid,nL,nMid);
    eR = EigsInInterval(T,Mid,R,nMid,nR);
    e = sort([eL;eR]);
end
    
    


  function lambda = Sturm(T,k,y,z)
% T is an nxn unreduced symmetric tridiagonal matrix and 1<=k<=n.
% lambda is the kth smallest eigenvalue of T.
n = length(T);
while abs(y-z)> eps*(abs(y)+abs(z))
    x = (y+z)/2;
    if EigsToLeft(T,x)>=k
        z = x;
    else
        y = x;
    end
end
lambda = (y+z)/2;

% -1 for many more calls to EigsToLeft than needed
    