   function [sigma] = bisect(d, b, i, tol)
%% function [sigma] = bisect(d, b, i, tol)
%% d   = diagonal of T
%% b   = co-diagonal
%% i   = compute i-th eigenvalue
%% tol = tolerance used for stopping
     b(1) = 0;
   n = length(d);
 %%-------------------- guershgorin
   tmin = d(n) - abs(b(n));
   tmax = d(n) + abs(b(n));
   for j=1:n-1
     rho = abs(b(j)) + abs(b(j+1));
     tmin = min(tmin, d(j)-rho);
     tmax = max(tmax, d(j)+rho);
   end
   tol = tol*(tmax-tmin);
   for iter=1:100
     sigma = 0.5*(tmin+tmax);
     count = sturm(d, b, sigma);
     if (count >= i) 
       tmin = sigma;
     else
       tmax = sigma;
     end
     if (tmax - tmin) < tol
       break
     end
   end
