function [mid_point,counter] = bisect2(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

%
counter=0;
x_guess_l=a;
x_guess_u=b;


while (abs(x_guess_u-x_guess_l)>tol && counter<imax)
    f_at_guess_l=feval(f,x_guess_l);
    f_at_guess_u=feval(f,x_guess_u);
    % Calculate mid-point
    mid_point=(x_guess_l+x_guess_u)/2;
    f_at_mid_point=feval(f,mid_point);

    % Output to screen
    fprintf('ITERATION %i \r',counter)
    fprintf('Lower guess = %10.6f \r',x_guess_l)
    fprintf('Upper guess = %10.6f \r',x_guess_u)
    fprintf('Mid-point = %10.6f \r \r',mid_point)
    
    % Decision 
    if f_at_mid_point*f_at_guess_u<=0
        x_guess_l=mid_point;
    else
        x_guess_u=mid_point;
    end
    counter=counter+1;
end