function [xr, n] = find_root_false_position(func, bracket, debug_flag)

%% ME 2016 Section A-B-C-D
%% Fall 2004
%%
%% AUTHOR:  Chris Paredis
%%
%% PURPOSE: find a root for the function defined in "func" inside the
%% bracket "bracket"
%%
%% ASSUMPTIONS:
%%         - The desired accuracy is expressed in absolute terms as 1e-6
%%         - A root with the desired accuracy will be found within 1000
%%         iterations; if not, an error is returned

%% we use the following variable names
%%   xu = upper limit of the bracket
%%   xl = lower limit of the bracket
%%   xr = next root approximation (inside bracket)
%%   fu = func(xu)
%%   fl = func(xl)
%%   fr = func(xr)

% desired accuracy of the algorithm
required_accuracy = 1e-6;
% exit criterion in case a root is not found quickly enough
max_iter = 3000;

% 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
% one function value must be positive, the other negative
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

% The main iteration loop of the algorithm
for n=3:max_iter
    % find the next root approximation
    xr = xu - fu*(xl-xu)/(fl-fu);
    fr = feval (func, xr);
    
    % check to see whether the debug_flag variable has been defined
    % (if more than 2 input variables exist, we should print debug statement)
    if nargin>2
        fprintf('xr = %g, fr = %g\n',xr,fr);
    end
    
    % determine the new bracket
    if fr == 0  % we could be lucky --> return immediately
        return;
    elseif 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
    
    % check whether the desired accuracy has been reached
    % this is the normal exit point of the function
    % NOTE: testing whether the function value is sufficiently small should
    % be used with care --> this is not safe if the function values are
    % small across a broad range of x-values
    if abs(fr) <= required_accuracy
        return;
    end
end

% if we reach this point, then we have iterated for too long
error('Maximum number of iteration (%d) exceeded\nCurrent root estimate: %g',max_iter, xr);