function [a, b, k, n] = min_golden_section(fun_expression, alpha, beta, epsilon, lambda)
% Golden section method for finding the local minimum of a function.
%
% fun_expr: (symbolic expression over x) The symbolic expression of the
%           objective function
% alpha:    (number) The starting point of the interval in which we seek
%           for minimum
% beta:     (number) The ending point of the interval in which we seek
%           for minimum
% epsilon:  (number) The epsilon value
%           **note:**
%             epsilon in not used in this method, but it is part of the
%             method calling interface.
% lambda:   (number) The lambda value (accuracy)
%
% return:
%  a:       (vector) Starting points of the interval for each iteration
%  b:       (vector) Ending points of the interval for each iteration
%  k:       (number) The number of iterations
%  n:       (number) The calls of objective function fun_expr
%

% Error checking
if alpha > beta || lambda <= 0
    error ('Input criteria not met')
end

% Init variables
gamma = 0.618;
a = alpha;
b = beta;
fun = matlabFunction(fun_expression);

% calculate x1, x2 of the first iteration, since the following iteration
% will not require to calculate both
k=1;
x_1 = a(k) + (1 - gamma)*(b(k) - a(k));
x_2 = a(k) + gamma*(b(k) - a(k));

f1 = fun(x_1);
f2 = fun(x_2);
while b(k) - a(k) > lambda    
    % set new search interval
    k = k + 1;
    if f1 < f2
        a(k) = a(k-1);
        b(k) = x_2;
        x_2 = x_1;
        f2 = f1;
        x_1 = a(k) + (1 - gamma)*(b(k) - a(k));
        f1 = fun(x_1);
    else
        a(k) = x_1;
        b(k) = b(k-1);
        x_1 = x_2;
        f1 = f2;
        x_2 = a(k) + gamma*(b(k) - a(k));
        f2 = fun(x_2);
    end
end

% Set objective function calls
n = k+1;