+#!/usr/bin/octave --silent
+
+function root = bisect(f, a, b, epsilon)
+ ## Find a root of the function `f` on the closed interval [a,b].
+ ##
+ ## It is assumed that there are an odd number of roots within [a,b].
+ ## The intermediate value theorem is used to determine whether or not
+ ## there is a root on a given interval. So, for example, cos(-2) and
+ ## cos(2) are both less than zero; the I.V.T. would not conclude that
+ ## cos(x) has a root on -2<=x<=2.
+ ##
+ ## If `f` has more than one root on [a,b], the smallest root will be
+ ## returned. This is an implementation detail.
+
+ ## Store these so we don't have to recompute them.
+ fa = f(a);
+ fb = f(b);
+
+ ## We first check for roots at a,b before bothering to bisect the
+ ## interval.
+ if (fa == 0)
+ root = a;
+ return;
+ end
+
+ if (fb == 0)
+ root = b;
+ return;
+ end
+
+ ## Bisect the interval.
+ c = (a+b)/2;
+
+ ## If we're within the prescribed tolerance, we're done.
+ if (b-c < epsilon)
+ root = c;
+ return;
+ end
+
+ if (has_root(fa,f(c)))
+ root = bisect(f,a,c,epsilon);
+ else
+ root = bisect(f,c,b,epsilon);
+ end
+end