function [p,yp,dp,dy,A,B,C,D] = golden(f,a,b,delta,epsilon)
% [p,yp,dp,dy,A,B,C,D] = golden(f,a,b,delta,epsilon)
% [p,yp,dp,dy] = golden(f,a,b,delta,epsilon)
% Golden ratio search for a minimum.
% f is the function, input.
% a is the left endpoint, input.
% b is the right endpoint, input.
% delta is the tolerance for the abscissas, input.
% epsilon is the tolerance for the ordinates, input.
% p is the abscissa of the minimum, output.
% dp is the ordinate of the minimum, output.
% dp is the error bound for p, output.
% dy is the error bound for yp, output.
% A is the vector of left endpoints, output.
% B is the vector of right endpoints, output.
% C is the vector of central values, output.
% D is the vector of central values, output.
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