code: def _densMoments_approx_higherorder_gaussxpolyInts(self,ll,ul,maxj):
          
          gaussxpolyInt= numpy.zeros((maxj,len(ul)))
          gaussxpolyInt[-1]= 1./numpy.sqrt(numpy.pi)\
              *(numpy.exp(-ll**2.)-numpy.exp(-ul**2.))
          gaussxpolyInt[-2]= 1./numpy.sqrt(numpy.pi)\
              *(numpy.exp(-ll**2.)*ll-numpy.exp(-ul**2.)*ul)\
              +0.5*(special.erf(ul)-special.erf(ll))
          for jj in range(maxj-2):
              gaussxpolyInt[-jj-3]= 1./numpy.sqrt(numpy.pi)\
              *(numpy.exp(-ll**2.)*ll**(jj+2)-numpy.exp(-ul**2.)*ul**(jj+2))\
              +0.5*(jj+2)*gaussxpolyInt[-jj-1]
          return gaussxpolyInt
docs: Calculate all of the polynomial x Gaussian integrals occuring 
          in the higher-order terms, recursively