code: def get_diff_coeff(hvec, n=1):
      
      hvec = np.array(hvec, dtype=np.float)
      acc = len(hvec)
      exp = np.column_stack([np.arange(acc)]*acc)
      a = np.vstack([hvec] * acc) ** exp
      b = np.zeros(acc)
      b[n] = factorial(n)
      return np.linalg.solve(a, b)
docs: Helper function to find difference coefficients of an
      derivative on an arbitrary mesh.
  
      Args:
          hvec (1D array-like): sampling stencil
          n (int): degree of derivative to find