function [x, N] = fibonacci_search( f, a, b, tol, maxfib, verbose )
  %---------------------------------------------------------------------------
  % Jonathan R. Senning <jonathan.senning@gordon.edu>
  % Gordon College
  % Written May 3, 1999
  % Revised December 5, 2000
  %
  % Implements the Fibonacci search algorithm as described on pages 552-555
  % of "Numerical Mathematics and Computing", 4th Edition, by Cheney & Kincaid,
  % Brooks/Cole, 1999.
  %---------------------------------------------------------------------------

  if ( nargin <= 5 )
    verbose = 0;
  end

  % Find the smallest Fibonacci number larger than 0.5 * (b - a) / tol.

  fib = [1, 1];
  for k = 3 : maxfib
    fib( k ) = fib( k - 1 ) + fib( k - 2 );
    if ( ( 0.5 * ( b - a ) / tol ) < fib( k ) )
      break;
    end
  end

  if ( ( 0.5 * ( b - a ) / tol ) < fib( k ) )
    N = k;
  else
    fprintf( 'Need more than %d Fibonacci numbers\n', maxfib );
    exit
  end

  if ( verbose ~= 0 )
    fprintf( 'N = %d; fib(%d) = %d\n', N, N, fib( N ) );
  end

  % Compute known error bound.  Here "d" plays the role of delta in the
  % algorithm on pages 553-554.

  d = 0.5 * ( b - a ) / fib( N );

  % Carry out the iterations necessary to compute our estimate of the minimum.
  % Here "D" plays the role of Delta in the algorithm on pages 553-554.
  % Similarly, "ap" and "bp" stand for a' and b' respectively.

  for k = N : -1 : 3
    D = ( b - a ) * fib( k - 2 ) / fib( k );
    ap = a + D;
    bp = b - D;
    if ( f( ap ) >= f( bp ) )
      a = ap;
    else
      b = bp;
    end
  end

  ap = 0.5 * ( a + b ) - 2 * d;
  bp = 0.5 * ( a + b ) + 2 * d;

  if ( f( ap ) >= f( bp ) )
    a = ap;
  else
    b = bp;
  end

  % Compute our final estimate of the solution.

  x = 0.5 * ( a + b );

end
