function [c,d,isroot] =bisect2(phi,c0,d0,tol)
%
% function [c,d,isroot] =bisect2(phi,c0,d0,tol)
%   Finds a root of a function phi which phi(c0)*phi(d0) <= 0.
%     phi--- String with name of MATLAB function
%     [c0,d0] -- Interval in which root of phi is believed to be
%     tol -- Tolerence at which root is to be calculated
%
%     Output parameters
%      [c,d] -- Interval of size about tol with a root of phi
%      isroot -- 0 if phi(c0)*phi(d0) > 0 else is equal to 1.
%
true = 1;
false =0;
c= c0;
d = d0;
x= ones(50);
phic = eval([phi,'(c)'])
phid = eval([phi,'(d)'])
if phic*phid < 0      
  count = 1; notdone=1;
  while notdone
    if count < 50
      mid = (c+d)/2
      phimid = eval([phi,'(mid)'])
      if phimid*phid < 0 
        c = mid;
        phic = phimid;
      else
        d = mid;
        phid = phimid;
      end
      x(count)=mid;
     count = count +1
     notdone = d-c > tol & abs(phimid) > tol;
    else 
      notdone=0;
    end
   end
   isroot = true;
   if abs(phimid) == 0
     c=mid;
     d=mid;
   end;
   x(1:count-1)
  else
   isroot = false ;
 end