function root = bisect(f_x,a,b,imax,tol)
% f       function to be evaluated 
% a,b     interval limits 
% imax    iteration maximum 
% tol     convergence criterion 
% root    position of the root 

% initialise parameters
x1=a
x3=b
f1=f_x(a)
f3=f_x(b)

% CHECK: verify that there is indeed a root in the interval [a b].
if f1*f3>0     
    error('No root in the specified interval') 
end

% Begin iteration loop i.e. for ndx =1:imax or while(ndx <imax)
counter=1
while abs(x3-x1)>tol && counter<=imax
% Find which half interval contains root – i.e. if-else construction 
x2=(x3+x1)/2;
f2=f_x(x2);
fprintf('midpoint = %10.6f\r',x2)

if f1*f2<0 
    % Root is in left half, so        
    d=(x2+x1)/2;     
    f3=f2;     
    x3=x2; 
else 
    % Root is in right half, so        
    d=(x2+x3)/2;     
    f1=f2;     
    x1=x2;
end
counter=counter+1
end
root=x2
end