function[xmin, fmin] = golden(xlow, xhigh, f)

maxit = 100;
es = 0.0001;

R = (sqrt(5) - 1)/2;

x1 = R*xlow + (1-R)*xhigh;
x2 = (1-R)*xlow + R*xhigh;

f1 = feval(f, x1);
f2 = feval(f, x2);

while (es < abs(xhigh - xlow))
    for i = 1:maxit-2
        if f1 < f2
            xhigh = x2;
            x2 = x1;
            f2 = f1;
            x1 = R*xlow + (1-R)*xhigh;
            f1 = feval(f, x1);
        else
            xlow = x1;
            x1 = x2;
            f1 = f2;
            x2 = (1-R)*xlow + R*xhigh;
            f2 = feval(f, x2);
        end
    end

    if f1 < f2
        xmin = x1;
        fmin = f(x1);
    else
        xmin = x2;
        fmin = f(x2);
    end
end



  