code: def derivative(f, t):
      
      dfdt = np.empty_like(f)
      if (f.ndim == 1):
          _derivative(f, t, dfdt)
      elif (f.ndim == 2):
          _derivative_2d(f, t, dfdt)
      elif (f.ndim == 3):
          _derivative_3d(f, t, dfdt)
      else:
          raise NotImplementedError("Taking derivatives of {0}-dimensional arrays is not yet implemented".format(f.ndim))
      return dfdt
docs: Fourth-order finite-differencing with non-uniform time steps
  
      The formula for this finite difference comes from Eq. (A 5b) of "Derivative formulas and errors for non-uniformly
      spaced points" by M. K. Bowen and Ronald Smith.  As explained in their Eqs. (B 9b) and (B 10b), this is a
      fourth-order formula -- though that's a squishy concept with non-uniform time steps.
  
      TODO: If there are fewer than five points, the function should revert to simpler (lower-order) formulas.