function [a, b, k] = bisection(fun, alpha, beta, epsilon, lambda)
%
% Detailed explanation goes here
%
%

% Error checking
if 2*epsilon >= lambda
    error ('Convergence criteria not met')
end

% Init output vectors
a = alpha;
b = beta;

k=1;

while b(k) - a(k) > lambda
    % bisect [a,b]
    mid = (a(k) + b(k)) / 2;
    x_1 = mid - epsilon;
    x_2 = mid + epsilon;
    
    % set new search reange
    k = k + 1;
    if fun(x_1) < fun(x_2)
        a(k) = a(k-1);
        b(k) = x_2;
    else
        a(k) = x_1;
        b(k) = b(k-1);
    end
end