function c = bisection(f,a,b,delta,epsilon,N)
%
% c = bisection(f,a,b,delta,epsilon,N)
%
% This is the bisection method applied to the inline function f on the
% interval [a,b].  The process will stop if either N iterations have
% completed or if the width of the "x" interval is less than 2*delta or
% if |f(c)| is less than epsilon where c is the midpoint of the
% iteratively refined interval [a,b].
%

% Troy Henderson
% Assistant Professor
% University of Mobile
% Updated 15 Jan 08

fa = f(a);

if (fa*f(b)>0)
	disp("\nIntermediate value theorem cannot be applied");
	disp("since f(a) and f(b) have the same sign.\n");
	error("Make sure f changes sign in [a,b].\n");
end

more off
for n=1:N
	d = (b-a)/2;
	if (d < delta)
		break;
	end
	c = a + d;
	fc = f(c);
%	printf("n = %d\t\tc = % .20f\tf(c) = % .8g\n",n,c,fc);
	printf("(%.6g,%.6g,%.6g)\n",a,c,b);
	if (abs(fc) < epsilon)
		break;
	end
	if (fa*fc > 0)
		a = c;
		fa = fc;
	else
		b = c;
	end
end

endfunction
