function root = bisect(f, a, b, imax, tol)
    % Input: 
    % f    - function to evaluate
    % a, b - interval limits
    % imax - maximum iterations
    % tol  - tolerance for convergence
    % Output:
    % root - position of the root
    
    x1 = a; 
    x3 = b;
    f1 = feval(f, x1);
    f3 = feval(f, x3);
    
    if f1 * f3 > 0
        error('No root in the specified interval')
    end
    
    for ndx = 1:imax
        x2 = (x1 + x3) / 2;
        f2 = feval(f, x2);
        
        if abs(f2) < tol
            root = x2;
            return;
        elseif f1 * f2 < 0
            x3 = x2;
            f3 = f2;
        else
            x1 = x2;
            f1 = f2;
        end
    end
    
    error('Maximum iterations reached without convergence')
end