--- /dev/null
+#!/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
--- /dev/null
+#!/usr/bin/octave --silent
+
+function has_root = has_root(fa, fb)
+ ## Use the intermediate value theorem to determine whether or not some
+ ## function has an odd number of roots on an interval. If the function
+ ## in question has an even number of roots, the result will be
+ ## incorrect.
+ ##
+ ## Call the function whose roots we're concerned with 'f'. The two
+ ## parameters `fa` and `fb` should correspond to f(a) and f(b).
+ ##
+
+ ## If either f(a) or f(b) is zero, the product of their signs will be
+ ## zero and either a or b is a root. If the product of their signs is
+ ## negative, then f(a) and f(b) are non-zero and have opposite sign,
+ ## so there must be a root on (a,b). The only case we don't want is
+ ## when f(a) and f(b) have the same sign; in this case, the product of
+ ## their signs would be one.
+ if (sign(fa) * sign(fb) != 1)
+ has_root = true;
+ else
+ has_root = false;
+ end
+end