class MoreMath::NewtonBisection

def solve(range = nil, n = 1 << 16, epsilon = 1E-16)
  if range
    x1, x2 = range.first.to_f, range.last.to_f
    x1 >= x2 and raise ArgumentError, "bad initial range #{range}"
  elsif range = bracket
    x1, x2 = range.first, range.last
  else
    raise DivergentException, "bracket could not be determined"
  end
  f = @function[x1]
  fmid = @function[x2]
  f * fmid >= 0 and raise DivergentException, "root must be bracketed in #{range}"
  root = if f < 0
           dx = x2 - x1
           x1
         else
           dx = x1 - x2
           x2
         end
  n.times do
    fmid = @function[xmid = root + (dx *= 0.5)]
    fmid < 0 and root = xmid
    dx.abs < epsilon or fmid == 0 and return root
  end
  raise DivergentException, "too many iterations (#{n})"
end