function [x, k] = golden_section( f, a, b, tol, maxiter, verbose )
  %---------------------------------------------------------------------------
  % Jonathan R. Senning <jonathan.senning@gordon.edu>
  % Gordon College
  % Written May 3, 1999
  % Revised May 4, 2005
  % Revised December 18, 2008 - Works with both Matlab and Octave
  %
  % Implements the golden section search as described on pages 555-556 of
  % "Numerical Mathematics and Computing", 4th Edition, by Cheney & Kincaid,
  % Brooks/Cole, 1999 except that y is calculated as an offset from b rather
  % than from a.
  %---------------------------------------------------------------------------

  if ( nargin <= 5 )
    verbose = 0;
  end
  
  % Set things up to begin the iteration
  r  = ( sqrt( 5 ) - 1 ) / 2;   % Golden ratio
  x = a + r * ( b - a );
  y = b - r * ( b - a );
  u = f( x );
  v = f( y );

  % Main loop
  k = 0;
  if ( verbose ~= 0 )
    fprintf( '%4d: Interval is [%f,%f]; y = %f, x = %f\n', k, a, b, y, x );
  end  
  while ( abs( b - a ) >= 2 * tol && k < maxiter )
    k = k + 1;
    if ( u > v )        % minimum must be in [a,x]
      b = x;
      x = y;
      u = v;
      y = b - r * ( b - a );
      v = f( y );
    else                % minimum must be in [y,b]
      a = y;
      y = x;
      v = u;
      x = a + r * ( b - a );
      u = f( x );
    end
    if ( verbose ~= 0 )
      fprintf( '%4d: Interval is [%f,%f]; y = %f, x = %f\n', k, a, b, y, x );
    end
  end

  % All done; best estimate of solution is midpoint of final interval
  x = 0.5 * ( a + b );

end
