function [xvect,xdif,fx,nit]=bisection(a,b,toll,nmax,fun,varargin)

% [xvect,xdif,fx,nit]=bisection(a,b,toll,nmax,fun,varargin)
%
% bisection method
%
%  >> INPUT <<
%
%  a, b       two values such that fun(a)*fun(b)<0       
%  toll       tolerance required
%  nmax       maximum number of iterations
%  fun        function f
%  varargin   additional variables of function fun
%             Warning, roots of function fun with respect to first
%             variable of function fun
%  
%
%  >> OUTPUT <<
%
%  xvect      vector containing all the iterates, the root is
%             xvect(end)
%  xdif       difference between to iterates
%  fx         values of fun in xvect
%  nit        number of iterations required to reach the 
%             prescribed tolerance
%  >> EXAMPLE <<
%
%  [xvect,xdif,fx,nit]=bisection(-1,0.9,1e-6,150,@sin);
%  [xvect,xdif,fx,nit]=bisection(-1,0.9,1e-6,150,'sin');
%  [xvect,xdif,fx,nit]=bisection(-1,0.9,1e-6,150,'atan2',1);
%  [xvect,xdif,fx,nit]=bisection(-1,0.9,1e-6,150,'F')  provided a function F.m
%  [xvect,xdif,fx,nit]=bisection(-1,0.9,1e-6,150,'F',a,b)  provided a
%  function F.m, , where F has three variables x, a and b.
%


sfa=sign(feval(fun,a,varargin{:}));
sfb=sign(feval(fun,b,varargin{:}));
if sfa*sfb==1
    error('la fonction a le même signe en a et b');
end
if sfa*sfb==0
    nit=0;
    xdif=[];
    fx=0;
    if sfa==0
        xvect=a;
    else
        xvect=b;
    end
else
    xvect=zeros(1,nmax);
    fx=xvect;
    xdif=xvect;
    fc=1;
    err=toll+1;
    nit=0;
    while (nit < nmax & err > toll & fc ~=0)
        nit=nit+1;
        c=(a+b)/2;
        fc=feval(fun,c,varargin{:});
        sfc=sign(fc);
        xvect(nit)=c;
        fx(nit)=fc;
        if (sfc*sfa>0)
            a=c;
            sfa=sfc;
        else
            b=c;
            sfa=-sfc;
        end;
        err=abs(b-a);
        xdif(nit)=err;
    end
    xvect=xvect(1:nit);
    fx=fx(1:nit);
    xdif=xdif(1:nit);
end