function root = bisect(fun,x1,x2)
% uses the bisection method to estimate a root
% of fun(x)=0 between x1 and x2.
% the method is iterated 20 times to provide
% 1 part in 2^20 of the original x1-x2 interval.
if fun(x1)*fun(x2) > 0
    root = 'initial estimates do not bracket the solution';
else
    for i=1:20
        xm = (x1+x2)/2;
        if fun(xm)*fun(x1) > 0
            x1 = xm;
        else
            x2 = xm;
        end
    end
    root = xm;
end

