function c = frobust(f,a,b)
% FROBUST  Putatively robust combination of secant method
%          and bisection.  To solve f(x) = 0.  Requires bracket
%          to get started.  Goal is to produce bracket with
%          length less than 10^-12.
% Easy example; much faster than bisection:
%   >> frobust(@(x) x-cos(x), 0.7, 0.8)
% Easy example but motivates subtle cases below:
%   >> frobust(@(x) abs(x) - 2, -1, 3)
% Hard example; complex behavior but success; plot this fcn to see:
%   >> frobust(@(x) abs(sin(4*x))- 0.2*x - 0.1,-1,3)
% Impossible example; THINKS it has a root, but it doesn't; evaluate f to see:
%   >> frobust(@(x) sin(1e16*x), -1, 1.05)

fa = feval(f,a);
fb = feval(f,b);
if sign(fa) == sign(fb), error('not a bracket!'), end
printf('                    initial bracket [%15.12f, %15.12f]\n', a, b)

while abs(b-a)/max([abs(a) abs(b)]) > 1e-12
  c = b - fb * (b - a) / (fb - fa)
  fc = feval(f,c);
  if (c >= a) & (abs(c-a) < 1e-12)
    % c is nearly same as a; try to move b very close to a
    [b, fb] = tryreduce_moveb(f,a,b,fb,min(a + 3e-12,b));
  elseif (c <= b) & (abs(b-c) < 1e-12)
    % c is nearly same as b; try to move a very close to b
    [a, fa] = tryreduce_movea(f,a,b,fa,max(b - 3e-12,a));
  end
  if (c > a) & (c < b) & (abs(fc) < min(abs(fa),abs(fb)))
    % accept c because it is in (a,b) and has better fcn value
    if (b - c) < 0.1 * (b - a)
      % if c is close to b, try to shrink the bracket by moving a
      [a, fa] = tryreduce_movea(f,a,b,fa,b - 3 * (b - c));
    elseif (c - a) < 0.1 * (b - a)
      % if c is close to a, try to shrink the bracket by moving b
      [b, fb] = tryreduce_moveb(f,a,b,fb,a + 3 * (c - a));
    end
    % now update bracket
    if sign(fc) == sign(fa)
      a = c;   fa = fc;
      printf('  secant step (move a); new bracket [%15.12f, %15.12f]\n', a, b)
    else
      b = c;   fb = fc;
      printf('  secant step (move b); new bracket [%15.12f, %15.12f]\n', a, b)
    end
  else
    % reject c; do bisection
    c = (a + b)/2;
    fc = feval(f,c);
    if sign(fc) == sign(fa)
      a = c;   fa = fc;
    else
      b = c;   fb = fc;
    end
    printf('        bisection step; new bracket [%15.12f, %15.12f]\n', a, b)
  end
end

    function [a, fa] = tryreduce_movea(f, a, b, fa, atry);
    fatry = feval(f,atry);
    if (sign(fatry) == sign(fa)) & (abs(fatry) < abs(fa))
      a = atry;   fa = fatry;
      printf('  reduce step (move a); new bracket [%15.12f, %15.12f]\n', a, b)
    end

    function [b, fb] = tryreduce_moveb(f, a, b, fb, btry);
    fbtry = feval(f,btry);
    if (sign(fbtry) == sign(fb)) & (abs(fbtry) < abs(fb))
      b = btry;   fb = fbtry;
      printf('  reduce step (move b); new bracket [%15.12f, %15.12f]\n', a, b)
    end
