docs: def derivative(func, x0, dx=1.0, n=1, args=(), order=3):
      '''Reimplementation of SciPy's derivative function, with more cached
      coefficients and without using numpy. If new coefficients not cached are
      needed, they are only calculated once and are remembered.
      '''
      if order < n + 1:
          raise ValueError
      if order % 2 == 0:
          raise ValueError
      weights = central_diff_weights(order, n)
      tot = 0.0
      ho = order >> 1
      for k in range(order):
          tot += weights[k]*func(x0 + (k - ho)*dx, *args)
      return tot/product([dx]*n)