In [2]:     mu1, mu2 = symbols('mu1 mu2', real=True, finite=True, bounded=True)
   ...:     sigma1, sigma2 = symbols('sigma1 sigma2', real=True, finite=True,
   ...:                                               bounded=True, positive=True)
   ...:     rate = Symbol('lambda', real=True, positive=True, bounded=True)
   ...:     def normal(x, mu, sigma):
   ...:             return 1/sqrt(2*pi*sigma**2)*exp(-(x-mu)**2/2/sigma**2)
   ...: 

In [3]:     def exponential(x, rate):
   ...:             return rate*exp(-rate*x)
   ...: 

In [4]: integrate(normal(x, mu1, sigma1), (x, -oo, oo), meijerg=True)
Out[4]: 1

In [5]: integrate(x*normal(x, mu1, sigma1), (x, -oo, oo), meijerg=True)
Out[5]: μ₁

In [6]: integrate(x**2*normal(x, mu1, sigma1), (x, -oo, oo), meijerg=True)
Out[6]: 
  2     2
μ₁  + σ₁ 

In [7]: integrate(normal(x, mu1, sigma1)*normal(y, mu2, sigma2),
   ...:                      (x, -oo, oo), (y, -oo, oo), meijerg=True)
Out[7]: 1

In [8]: integrate(x*y*normal(x, mu1, sigma1)*normal(y, mu2, sigma2),
   ...:                      (x, -oo, oo), (y, -oo, oo), meijerg=True)
Out[8]: μ₁⋅μ₂

In [9]: integrate(x**2*exponential(x, rate), (x, 0, oo), meijerg=True)
Out[9]: 
2 
──
 2
λ 

In [12]: k = Symbol('k', integer=True, positive=True)

In [13]: chisquared = 2**(-k/2)/gamma(k/2)*x**(k/2-1)*exp(-x/2)

In [14]: powsimp(integrate(chisquared, (x, 0, oo), meijerg=True))
Out[14]: 1

In [16]: simplify(integrate(x**2*chisquared, (x, 0, oo), meijerg=True))
Out[16]: k⋅(k + 2)