function root = bisect1(fun,x1,x2,varargin)
% uses the bisection method to estimate a root
% of fun(x)=0 between x1 and x2.
% number of iterations is 20 by default
% 1 part in 2^20 of the original x1-x2 interval, but
% additional argument can specify iterations.
nargin = length(varargin);
if nargin == 0  
    maxit = 20;  % here if no extra argument
else
    maxit = cell2mat(varargin(1));  % here with
end                                 % extra argument
if fun(x1)*fun(x2) > 0
    root = 'initial estimates do not bracket the solution';
else
    for i=1:maxit
        xm = (x1+x2)/2;
        if fun(xm)*fun(x1) > 0
            x1 = xm;
        else
            x2 = xm;
        end
    end
    root = xm;
end

