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

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

% Init output vectors
a = alpha;
b = beta;
dfun = matlabFunction(diff(fun_expression));

k=1;
while b(k) - a(k) > lambda
    % bisect [a,b]
    x_mid = (a(k) + b(k)) / 2;
    
    % set new search interval
    k = k + 1;
    df = dfun(x_mid);
    if df < 0
        a(k) = x_mid;
        b(k) = b(k-1);
    elseif df > 0
        a(k) = a(k-1);
        b(k) = x_mid;
    else % df == 0
        a(k) = x_mid;
        b(k) = x_mid;
        break;
    end
end