function interval = hybrid(f,a,b,delta)
% interval = hybrid(f,a,b,delta)
%
%   find interval [x,y] with f(x)*f(y)<=0
%   requires a,b with f(a)*f(b) < 0
%
% algorithm with following three properties:
%  * gives interval containing zero after each step ("bracketing")
%  * for continuous f, length of interval converges to zero with order 1
%  * for f with f'' continuous, length of interval converges to zero with
%    superlinear order (1.618).

% T. von Petersdorff 04/08/2004
% based on hybrid secant-bisection method from Novak, Ritter, Wozniakowski (Math. Comp. 1995)

% This is just a demonstration implementation of an algorithm for which one can prove 
% the three properties above. If e.g. delta is chosen too small with respect to the machine epsilon 
% it will not terminate.
%
% For practical applications use fzero (based on Dekker-Brent algorithm) which
% is more efficient and robust.

fa = feval(f,a); 
fb = feval(f,b);
if fa==0
  interval = [a,a];
elseif fb==0
  interval = [b,b];
elseif sign(fa)==sign(fb)
  error('f must have different signs at a,b')
end

x = NaN; fx = NaN; l1 = NaN; l2 = NaN;

k = 0;

while abs(b-a) > delta
  k = k + 1;
  if k <= 2                                  %  use regula falsi 
    type = 'regula falsi';
    xnew = b - fb*(b-a)/(fb-fa);
  else
    accepted = 0;
    if abs(b-a) <= l3/2 | k==3               %  length reduced by .5 after 3 steps or 1st step after r.f.
      type = 'secant';                       %  consider secant step
      d = - fx*(x-xo)/(fx-fxo);
      xnew = x + d;
      if xnew==x                             %  make sure xnew and x are different machine numbers
        xnew = x*(1+sign(d)*sign(x)*eps);
      end
      if sign(xnew-a)*sign(xnew-b) <  0      %  xnew is in (a,b)
        accepted = 1;
      end
    end
    if ~accepted                             %  if not accepted
      type = 'bisection';
      xnew = (a+b)/2;                        %  use bisection step
      k = 0;
    end
  end
  xo = x; x = xnew; fxo = fx;
  l3 = l2; l2 = l1; l1 = abs(b-a);           %  save old lengths
  fx = feval(f,x);                           %  continue with subinterval which has sign change
  if fx==0                                   %  use [x,x]
    a = x; b = x; fa = 0; fb = 0;
  elseif sign(fa)*sign(fx)<0                 %  use [a,x]
    b = x; fb = fx;
  else                                       %  use [x,b]
    a = x; fa = fx;
  end  
  fprintf('[ %-20.17g, %-20.17g ], L = %-10.5g, %s\n', a, b, abs(b-a), type);
end
interval = [a,b];
  
