function [alpha] = gsm(f,a,b,eps,N)

%  golden section search on the function f. 

% You should define the f as function myfun.

%  No more than N function evaluation are done.

%  When b-a < eps, the iteration stops.
%
%     [alpha] = gsm('myfun',a,b,epsilon,N)
%  Example :
%     [alpha] = gsm('myfun',0,3,0.01,20)
%

clc
L = b-a;
x1 = a + .382*L;
fx1 = feval(f,x1);
x2 = a + .618*L;
fx2 = feval(f,x2);
fprintf('--------------------------------------------------------\n');
fprintf('    x1         x2         f(x1)       f(x2)       b - a\n');
fprintf('--------------------------------------------------------\n');
fprintf('%.3e  %.3e  %.3e  %.3e  %.3e\n', x1, x2, fx1, fx2, b-a);
for i = 1:N-2
   if fx1 < fx2
      b = x2;
      x2 = x1;
      fx2 = fx1;
      L=b-a;
      x1 = a + .382*L;
      fx1 = feval(f,x1);
   else
      a = x1;
      x1 = x2;
      fx1 = fx2;
      L=b-a;
      x2 = a + .618*L;
      fx2 = feval(f,x2);
   end;
   fprintf('%.3e  %.3e  %.3e  %.3e  %.3e\n', x1, x2, fx1, fx2, b-a);
   if (abs(b-a) < eps)
      fprintf('succeeded after %d steps\n', i);
      alpha = (b+a)/2;
      fprintf('alpha = %.4e\n',alpha);
      return;
   end;
end;
fprintf('failed requirements after %d steps\n', N);


