% Matlab code bisect.m 
% For "Applied Numerical Linear Algebra", Algorithm 1.1 and Question 1.21
% James Demmel, Aug 1995, modified May 1997
%
% function x = bisect(coeff,a,b,tol)
%
% Find a zero of a polynomial p(x)=0 
% given an interval [a,b] where p(a)*p(b) < 0.
%
% Inputs:
%   coeff = vector of polynomial coefficients
%           p(x) = sum from i=0 to n of coeff(i+1)*x^(n+1-i)
%                  where n = length(coeff)+1
%   a,b   = interval containing (at least) one zero
%   tol   = stopping tolerance. stop bisecting if
%           the interval containing a zero has width < 2*tol
%
% Outputs:
%   x = [xlow,xhigh] = interval containing a zero
%  
function x = bisect(coeff,a,b,tol)
  xlow = a; 
  xhigh = b;
  plow  = polyval(coeff,xlow);
  phigh = polyval(coeff,xhigh);
  while xhigh - xlow > 2*tol,
    xmid = (xlow + xhigh)/2;
    pmid = polyval(coeff,xmid);
    if pmid*plow < 0,
       xhigh = xmid;
       phigh = pmid;
    elseif pmid*phigh < 0,
       xlow = xmid;
       plow = pmid;
    else
       xlow = xmid;
       xhigh = xmid;
    end
  end
  x = [xlow, xhigh];