function x = LambdaSmall(a,c)
% a is a column n-vector and c is a column (n-1)-vector with nonzero entries.
% The tridiagonal matrix T(a,c) has at least one positive and one negative eigenvalue.
% x is the eigenvalue of T(a,c) that is closest to zero.

n = length(a); 
normT = max(abs(a) + abs([c;0]) + abs([0;c]));

k = s(0,a,c);
% k = number of eigenvalues <=  0. Thus, if we order
% the eigenvalues lambda(1)< ... < lambda(n), then we see that 
% x = lambda(k) or lambda(k+1).

% Let's get an interval [L,0] with the property that s(L,a,c) = k-1
L = -normT; sL = 0;
R = 0; sR = k;
while sL < k-1
    m = (L+R)/2; sm = s(m,a,c);
    if sm == k
        R = m;
    else
        L = m; sL = sm;
    end
end
% We now have L <= lambda(k) <= R and no other eigs in [L,R],so use fzero...
lambda_k = fzero(@(x) pEval(x,a,c),[L R]);
fprintf('\n L = %20.15f\n x = %20.15f\n R = %20.15f\n',L,lambda_k,R)

% Let's get an interval [0,R] with the property that s(R,a,c) = k+1
L = 0; sL = k;
R = normT; sR = n;
while k+1<sR
    m = (L+R)/2; sm = s(m,a,c);
    if sm == k
        L  = m;
    else
        R = m; sR = sm;
    end
end
% We now have L <= lambda(k+1) <= R and no other eigs in [L,R] so use fzero...
lambda_kp1 = fzero(@(x) pEval(x,a,c),[L,R]);
fprintf('\n L = %20.15f\n x = %20.15f\n R = %20.15f\n',L,lambda_kp1,R)

% Return the smaller...
if abs(lambda_k)<abs(lambda_kp1)
    x = lambda_k;
else
    x = lambda_kp1;
end

    
function p = pEval(x,a,c)
% p = det(T(a,c)-xI) via the 3-term recurrence..
n = length(a); pkm2 = 1;pkm1 = a(1)-x;
for k=2:n
    pk = (a(k)-x)*pkm1 - c(k-1)^2*pkm2;
    pkm2 = pkm1; 
    pkm1 = pk;
end
p = pkm1;

function sigma = s(x,a,c)
% sigma is the number of sign changes in the sequence
% {p0,p1,...,pn} where pk = det(T(1:k,1:k)-xI) where T = T(a,c).

n = length(a); pkm2 = 1; pkm1 = a(1)-x;
if pkm1>0
    Sign = 1; sigma = 0;
else
    Sign = -1; sigma = 1;
end
% Sign is +1 if the current pk is positive and -1 if not.
% sigma keeps track of the number of sign changes so far.
for k=2:n
    pk = (a(k)-x)*pkm1 - c(k-1)^2*pkm2;
    pkm2 = pkm1; 
    pkm1 = pk;
    if (pk<0 && Sign==1) || (pk>0 && Sign==-1) || pk==0
        sigma = sigma+1; Sign = -Sign;
    end
end