function x = IntervalBisection(fun,x,tol,varargin)

% x = IntervalBisection(f,x,tol,varargin)
%
% applies the interval bisection method for root-finding
%
%   f           interval function
%   x           at input:  search interval
%               at output: interval enclosing the roots
%   tol         relative error 
%   varargin    additional arguments for f

while max(relerr(x))>tol
    x = subdivide1D(x);
    f = x; 
    for k=1:length(x)
        f(k) = feval(fun,x(k),varargin{:});
    end
    rem = not(in(zeros(size(f)),f));
    f(rem) = []; x(rem) = [];
end
x=infsup(min(inf(x)),max(sup(x)));

return