[S,E,G]=golden(f,a,b,delta,epsilon)

function [S,E,G]=golden(f,a,b,delta,epsilon)

%Input  - f is the object function input as a string 'f'
%       - a and b are the endpoints of the interval
%       - delta is the tolerance for the abscissas
%       - epsilon is the tolerance for the ordinates
%Output - S=(p,yp) contains the abscissa p and
%         the ordinate yp of the minimum
%       - E=(dp,dy) contains the error bounds for p and yp
%       - G is an n x 4 matrix: the kth row contains [ak ck dk bk];
%         the values of a, c, d, and b at the kth iteration

% NUMERICAL METHODS: MATLAB Programs
%(c) 1999 by John H. Mathews and Kurtis D. Fink
%To accompany the textbook:
%NUMERICAL METHODS Using MATLAB,
%by John H. Mathews and Kurtis D. Fink
%ISBN 0-13-270042-5, (c) 1999
%PRENTICE HALL, INC.
%Upper Saddle River, NJ 07458

r1=(sqrt(5)-1)/2;
r2=r1^2;
h=b-a;
ya=feval(f,a);
yb=feval(f,b);
c=a+r2*h;
d=a+r1*h;
yc=feval(f,c);
yd=feval(f,d);
k=1;
A(k)=a; B(k)=b; C(k)=c; D(k)=d;

while (abs(yb-ya)>epsilon)|(h>delta)
   k=k+1;
   if (yc<yd)
      b=d;
      yb=yd;
      d=c;
      yd=yc;
      h=b-a;
      c=a+r2*h;
      yc=feval(f,c);
   else
      a=c;
      ya=yc;
      c=d;
      yc=yd;
      h=b-a;
      d=a+r1*h;
      yd=feval(f,d);
   end
   A(k)=a; B(k)=b; C(k)=c; D(k)=d;
end

dp=abs(b-a);
dy=abs(yb-ya);
p=a;
yp=ya;

if (yb<ya)
   p=b;
   yp=yb;
end

G=[A' C' D' B'];
S=[p yp];
E=[dp dy];