%
%   BISECTION.M  Routine to implement the Bisection Method
%
%   Inputs:  f   = function handle for user-supplied function.  
%                  The function f is of the form:  function y = f(x)
%            a,b = lower and upper limits (where b > a)
%            tol = convergence criterion  (tol > 100*eps)
%            M   = maximum number of iterations (M >= 2)
%        display = intermediate results are displayed if > 0                     .
%
%   Outputs:   x = estimated root of f(x) = 0 
%              k = number of iterations performed  (note that the number 
%                  of function evaluations is k + 2)
%   
%   Note that a and b must be selected such that f(a)*f(b) < 0 to ensure that  
%   there is a root in [a,b]
%
%   File prepared by J. R. White, UMass-Lowell (last update: Nov. 2017)
%
      function [x,k] = bisection(f,a,b,tol,M,display)
%
%   check bounds, tolerance, # iterations, and display switch
      fa = f(a);   fb = f(b);
      if fa*fb >= 0
        fprintf ('\n  Warning -- In bisection.m it is required that f(a)f(b) < 0.\n\n');
        return
      end
      if tol < 100*eps,  tol = 100*eps;  end
      if M < 2,  M = 2; end
      if display > 0
        fprintf('\n\n Intermediate edit from the Bisection Method:  \n\n');
        fprintf('    k       a         f(a)        b         f(b)       xr       f(xr)    sign(fa*fr)\n');
      end
%
%   find root 
      k = 1;      err = tol + 1;
      while (err > tol) && (k <= M)
        xr = a + (b-a)/2;   fr = f(xr);
        if display > 0
          fprintf('  %3i %10.6f %10.6f %10.6f %10.6f %10.6f %10.6f  %6i  \n', ...
                  k,a,fa,b,fb,xr,fr,sign(fa*fr));
        end
        if fa*fr < 0
          b = xr;   fb = fr;
        else 
          a = xr;   fa = fr;
        end
        err = abs(fr);   k = k + 1;	   
      end
      x = xr;  k = k-1;   % save outputs
      if k == M 
        fprintf('\n\n  Warning: The maximum number of iterations was reached!!! \n');
      end
%
%   end of function
