% Demonstration of Bisection Root Finding of a function

function y = fcn(x)
y = 9.8 * 68.1 * (1-exp(-10*x/68.1))/x - 40;
end


xl = 1;   % 'Enter lower bound of root bracket
xu = 4;   % 'Enter upper bound of root bracket
maxerror = 0.01;  % 'Enter maximum error 
maxit = 100;  % 'Enter maximum number of iterations



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

write("xl  ","xu  ","xr  ","error  ");

% main loop until interval width small enough

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

  count = count + 1;
  xr = (xl + xu)/2;
 
   if xr ~= 0                            % ~= not equal
    actual_error = abs((xu - xl)/(xu + xl)) * 100;
  end
  write(xl,"  ",xu,"  ",xr,"  ",actual_error);
  test = fcn(xl) * fcn(xr);  % form test product

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