function [xStar,fStar,nEvals,aF,bF] = GlobalNewton(f,derf,a,b,tolx,tolf,nEvalsMax)
% f handle to a function f(x).
% derf handle to its derivative f'(x).
% a,b [a,b] is a bracketing interval, i.e. f(a)*f(b)<=0.
% tolx,tolf Nonnegative termination criteria.
% nEvalsMax Maximum number of f-evaluations.
%
% xStar An approximate zero of f.
% fStar The value of f at xStar.
% nEvals The number of evaluations required.
% aF,bF [aF,bF] is the final bracketing interval.
%
% Generally xStar is within tolx of a true zero. However,the iteration
% terminates if |f(x)|<= tolf or if the number of f-evaluatios is about to
% exceed nEvalMax.
fa = f(a);
fb = f(b);
if fa*fb>0
disp('Initial interval not bracketing.')
return
end
xc = a;
fprintf('Initial: %20.15f\n',xc)
fc = f(xc);
dfc = derf(xc);
nEvals = 1;
while nEvals==1 || (abs(a-b) > tolx && abs(fc) > tolf && nEvals < nEvalsMax )
%[a,b] brackets a root and xc = a or xc = b.
if StepIsIn(xc,fc,dfc,a,b)
%Take Newton Step
xc = xc-fc/dfc;
fprintf('Newton: %20.15f\n',xc)
else
%Take a Bisection Step:
xc = (a+b)/2;
fprintf('Bisection: %20.15f\n',xc)
end
fc = f(xc);
dfc = derf(xc);
nEvals = nEvals+1;
% Update the bracketing interval...
if fa*fc<=0
% There is a root in [a,xc]. Bring in right endpoint.
b = xc;
fb = fc;
else
% There is a root in [xc,b]. Bring in left endpoint.
a = xc;
fa = fc;
end
end
aF = a;
bF = b;
xStar = xc;
fStar = fc;

function ok = StepIsIn(x,fx,fpx,a,b)
% ok = StepIsIn(x,fx,fpx,a,b)
% Yields 1 if the next Newton iterate is in [a,b] and 0 otherwise.
% x is the current iterate, fx is the value of f at x, and fpx is
% the value of f' at x.

if fpx > 0
ok = ((a-x)*fpx <= -fx) & (-fx <= (b-x)*fpx);
elseif fpx < 0
ok = ((a-x)*fpx >= -fx) & (-fx >= (b-x)*fpx);
else
ok = 0;
end