code: def differential(poly, diffvar):
      
      poly = Poly(poly)
      diffvar = Poly(diffvar)
  
      if not chaospy.poly.caller.is_decomposed(diffvar):
          sum(differential(poly, chaospy.poly.caller.decompose(diffvar)))
  
      if diffvar.shape:
          return Poly([differential(poly, pol) for pol in diffvar])
  
      if diffvar.dim > poly.dim:
          poly = chaospy.poly.dimension.setdim(poly, diffvar.dim)
      else:
          diffvar = chaospy.poly.dimension.setdim(diffvar, poly.dim)
  
      qkey = diffvar.keys[0]
  
      core = {}
      for key in poly.keys:
  
          newkey = np.array(key) - np.array(qkey)
  
          if np.any(newkey < 0):
              continue
          newkey = tuple(newkey)
          core[newkey] = poly.A[key] * np.prod(
              [fac(key[idx], exact=True) / fac(newkey[idx], exact=True)
               for idx in range(poly.dim)])
  
      return Poly(core, poly.dim, poly.shape, poly.dtype)
docs: Polynomial differential operator.
  
      Args:
          poly (Poly) : Polynomial to be differentiated.
          diffvar (Poly) : Polynomial to differentiate by. Must be decomposed. If
                  polynomial array, the output is the Jacobian matrix.
  
      Examples:
          >>> q0, q1 = chaospy.variable(2)
          >>> poly = chaospy.Poly([1, q0, q0*q1**2+1])
          >>> print(poly)
          [1, q0, q0q1^2+1]
          >>> print(differential(poly, q0))
          [0, 1, q1^2]
          >>> print(differential(poly, q1))
          [0, 0, 2q0q1]