function [a, b, k] = golden_section(fun, alpha, beta, lambda)
%


% Error checking
if lambda <= 0
    error ('Convergence criteria not met')
end

% Init variables
gamma = 0.618;
a = alpha;
b = beta;


% 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));

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