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
