function [x,X]=bisect(f,ab,tol)
% bisect: Find approx. of a zero of function f by bisection
% Inputs: f function handle, function whose zero is sought.
%         ab: Initial interval [a b]
%         tol: error tolerance, default: eps
% Result: x: Approksimation of a zero of f
%         X: The whole itaration sequence ( x is X(end)))
% Example:
%        [x,X] = bisect(@sin,[pi/2 3*pi/2]) 
%         fplot(@sin,[pi/2 3*pi/2]);grid on;shg
%         length(X)    % Gives nr. of iterations.
%
%         bisect(@sin,[2,2.5])  % Should cause error.

a=ab(1);b=ab(2);
if sign(f(a))==sign(f(b))
    error('f(a) and f(b) have same signs')
end    
if nargin==2 
    tol=eps;
end    
 X=zeros(1,55);  % Allocate X-vector for efficiency (52 + small)
                 % Note: eps=2^{-52}, nr. of steps is independent of f
 k=0;
 while((abs(b-a)) > (tol*abs(b)) & (abs(b-a))> tol)
     k=k+1;
     x=(a+b)/2; X(k)=x;  
     if sign(f(x)) == sign(f(b))
       b=x;
     else
       a=x;
     end
 end
   x=(a+b)/2;X=[X x]; 
