bisect.m

% Demonstration of Bisection Root Finding of a function
%
% The function is defined in the file 'fcn.m' as follows:
%      function y = fcn(x)
%      y = 9.8 * 68.1 * (1-exp(-10*x/68.1))/x - 40;

clear;   

      xl = input ('Enter lower bound of root bracket:  ');
      xu = input ('Enter upper bound of root bracket:  ');
maxerror = input ('Enter maximum error (percent):      ');
   maxit = input ('Enter maximum number of iterations: ');

fplot('fcn',[xl xu]); grid on; hold on;

count = 0;          % iteration counter
actual_error = 1;   % to force entry to while loop

fprintf('\n  #         xl         xu         xr      error\n\n');

% main loop until interval width small enough

while (actual_error > maxerror) & (count < maxit)

  count = count + 1;
  xr = (xl + xu)/2;
  plot(xr,fcn(xr),'ro');

  if xr ~= 0                            % ~= not equal
    actual_error = abs((xu - xl)/(xu + xl)) * 100;
  end

  fprintf('%3g %10g %10g %10g %10.4f\n',count,xl,xu,xr,actual_error)

  test = fcn(xl) * fcn(xr);  % form test product

  if test == 0
    actual_error = 0;
  elseif test < 0
    xu = xr;      % root is below xr
  else        % i.e. test > 0
    xl = xr;      % root is above xr
  end

end