function [b, k, vec] = dekker(f, a, b)
   %DEKKER dekker(f,a,b) reduces the interval [a,b] to a tiny
   % interal on which the function f(x) changes sign.
   % Returns one end point of that interval, the number
   % of function evaluations, and the vector of approximations

   vec = [];
   fa = f(a);
   fb = f(b);
   if sign(fa) == sign(fb)
      error('f(x) must change sign on the interval [a,b].')
   end
   k = 2;

   % a is the previous value of b and [b, c] always contains the zero.
   c = a; 
   fc = fa;
   
   while k < 50
      if sign(fb) == sign(fc)
         c = a; 
         fc = fa;
      end

      % Swap to insure f(b) is the smallest value so far.
      if abs(fc) < abs(fb)
         a = b; 
         fa = fb;
         b = c; 
         fb = fc;
         c = a; 
         fc = fa;
      end

      % Midpoint.
      m = (b + c)/2;
      if abs(m - b) <= eps(abs(b))
         return   % Exit from the loop and the function here.
      end

      % p/q is the the secant step.
      p = (b - a)*fb;
      if p >= 0
         q = fa - fb;
      else
         q = fb - fa;
         p = -p;
      end
      
      % Save this point.
      a = b; 
      fa = fb;

      % Choose next point.
      if p <= eps(q)
         % Minimal step.
         b = b + sign(c - b)*eps(b);
         fb = f(b);
      elseif p <= (m - b)*q
         % Secant.
         b = b + p/q;
      else
         % Bisection.
         b = m;
      end
      fb = f(b);
      k = k + 1;
      vec = [vec b];
   end
end

