function [xnew,its]=bisect(f,a,b,tol)

%BISECT   bisection method
%   [xnew,its]=bisect(f,a,b,tol) uses the bisection method for rootfinding. 
%   It returns a root, xnew, of f(x) in the interval [a,b] to within a 
%   tolerance tol. The number of iterations, its, is also returned.
%
% MMSC, October 2021
fa=f(a);
fb=f(b);
its=0;

if abs(fa)<tol
    % a is a root
    xnew=a;
elseif abs(fb)<tol
    % b is a root
    xnew=b;
elseif fa*fb>0
    error('f(a) and f(b) do not have opposite sign')
else
    % do bisection
    while (abs(b-a)>tol && its<2000)
        c=(a+b)/2;
        fc=f(c);
        if fc*fa>0
            % f(a) and f(c) have the same sign
            % so replace a with c
            a=c;
            fa=fc;
        else
            % replace b with c
            b=c;
            % fb=fc; % fb not actually needed
        end
        its=its+1;
    end
    xnew=c;
end