code: def differential(P, Q):
      
      P, Q = Poly(P), Poly(Q)
  
      if not chaospy.poly.is_decomposed(Q):
          differential(chaospy.poly.decompose(Q)).sum(0)
  
      if Q.shape:
          return Poly([differential(P, q) for q in Q])
  
      if Q.dim>P.dim:
          P = chaospy.poly.setdim(P, Q.dim)
      else:
          Q = chaospy.poly.setdim(Q, P.dim)
  
      qkey = Q.keys[0]
  
      A = {}
      for key in P.keys:
  
          newkey = numpy.array(key) - numpy.array(qkey)
  
          if numpy.any(newkey<0):
              continue
  
          A[tuple(newkey)] = P.A[key]*numpy.prod([fac(key[i], \
              exact=True)/fac(newkey[i], exact=True) \
              for i in range(P.dim)])
  
      return Poly(B, P.dim, P.shape, P.dtype)
docs: Polynomial differential operator.
  
      Args:
          P (Poly):
              Polynomial to be differentiated.
          Q (Poly):
              Polynomial to differentiate by. Must be decomposed. If polynomial
              array, the output is the Jacobian matrix.