function [root,counter] = bisect(f, a, b, imax, tol)
    % Input:
    %   f     - function to be evaluated
    %   a, b  - interval limits
    %   imax  - maximum number of iterations
    %   tol   - tolerance for convergence
    % Output:
    %   root  - position of the root

% Set iteration counter

    x1_guess = a;
    f1 = feval(f, x1_guess);
    x3_guess = b;
    f3 = feval(f, x3_guess);
    
    % [a, b]
    if f1 * f3 > 0
        error('No root in the specified interval')
    end
    
    % loop
    for ndx = 1:imax
        % 
        x2 = (x1_guess + x3_guess) / 2;
        f2 = feval(f, x2);

        % convergence
        if abs(f2) < tol || (x3_guess - x1_guess) / 2 < tol
            root = x2;
            return
        end

        % the root
        if f1 * f2 < 0
            x3_guess = x2; % root is in [x1, x2]
            f3 = f2;
        else
            x1_guess = x2; % root is in [x2, x3]
            f1 = f2;
        end
    end
    
    % If max iterations are reached
    root = (x1_guess + x3_guess) / 2;
    warning('Maximum iterations reached, result may not be within tolerance')
end