function [xr, n] = find_root_bisection(func, bracket, rel_acc)

%FIND_ROOT_BISECTION     finds a root using the bisection algorithm
%   [XR, N] = FIND_ROOT_BISECTION( FUNC, BRACKET, ACCURACY)
%
%   Inputs:
%      FUNC:       function handle for the function whose root is to be found
%      BRACKET:    a vector of length 2 indicating the initial bracket
%      ACCURACY:   the desired relative accuracy
%
%   Outputs:
%      XR:         The estimate of the root
%      N:          The number of function evaluations needed to determine
%                  the root  

%% ME 2016
%% Fall 2006
%% AUTHOR:  Chris Paredis
%%

% check the size of the bracket
if numel(bracket) ~= 2
    error('Initial bracket must be a vector of length 2');
end

% initialize the bracket variables
xu = bracket(2);
xl = bracket(1);
fu = feval (func, xu);
fl = feval (func, xl);
n = 2;

% check whether the bracket is valid
if fu == 0
    xr = xu;
    return;
elseif fl == 0;
    xr = xl;
    return;
elseif fu*fl >= 0
    error('You must provide an interval that brackets the root');
end

% exit criterion in case a root is not found quickly enough
max_func_eval = 3000;

% The main iteration loop of the algorithm
for n=3:max_func_eval
    % find the next root approximation
    xr = (xl + xu)/2 ;
    fr = feval (func, xr);
    
    % check whether the desired accuracy has been reached
    % this is the normal exit point of the function
    % NOTE: xrnew - xrold is the same as (xu-xl)/2
    if (fr == 0) || (abs(xu-xl)/2 <= rel_acc*xr)
        return;
    elseif (xl == xr) || (xu == xr)
        % NOTE: it is important to check whether xr is equal to one of the
        % bracket bounds --- otherwise the algorithm could get stuck with
        % the same bracket without ever reaching the required accuracy
        % EXAMPLE:
        % xl = 1;
        % xu = 1 + eps;
        % xr = (xl + xu) / 2; -->  oops: results in 1 (same as xl)        
        warning('The desired accuracy cannot be obtained');
        return;
    end
    
    % determine the new bracket
    if fr*fu < 0  
        % if sign of fr is opposite of fu, xr becomes the lower limit
        xl = xr;
        fl = fr;
    else
        % if sign of fr is opposite of fl, xr becomes the upper limit
        xu = xr;
        fu = fr;
    end
end

% if we reach this point, then we have iterated for too long -- however,
% this should NEVER happen
error('Too many iterations -- NOTE: this should never happen');