function root = bisect(f,a,b,imax,tol)
% This Function uses the bisection method to root solve equation 'f'...
% bewteen 'a' and 'b' to a tolerence of 'tol' with max iterations 'imax'.
% Input: f function to be evaluated
% a,b interval limits
% imax iteration maximum
% tol convergence criterion
% Output: root position of the root
% initialise parameters, i.e.

x1=a;           % x1 is lower guess for x
x3=b;           % x3 is upper guess for x
% x2 is midpoint for x
f1=feval(f,x1); 
f3=feval(f,x3);
if f1*f3>0
% CHECK: verify that there is indeed a root in the interval [x1(a) x3(b)].
error('No root in the specified interval')
end

diff = inf;
counter = 0;
while diff>tol && counter<imax
   f1 = feval(f,x1);
   f3 = feval(f,x3);
    
    x2 = (x3+x1)/2;
    f2 = feval(f,x2);

    % fprintf('ITERATION %i \r',counter)
    % fprintf('Lower guess = %.7f \r',x1)
    % fprintf('Upper guess = %.7f \r',x3)
    % fprintf('Mid-point = %.7f \r \r',x2)

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

    root = x2;
    fprintf('ITERATION %i \r',counter)
    fprintf('Root Estimate (7d.p) = %.7f \r\r',root)
end