The integral of the product of two normal distribution densities can be exactly solved, as shown here for example.
I'm interested in compute the following integral (for a generic $n \in \mathbb{N}$): $$ I_n = \int_{-\infty}^{\infty} \frac{1}{\sqrt{2 \pi} \sigma} e^{-\frac{x^2}{2 \sigma^2}} \frac{1}{\sqrt{2 \pi} \rho} e^{-\frac{(x-\mu)^2}{2 \rho^2}} x^n\mathrm{d}x $$
Is it possible to extend this result and solve $I_n$ for a generic $n \in \mathbb{N}$?
$\def\m{\mu} \def\p{\pi} \def\s{\sigma} \def\f{\varphi} \def\r{\rho} \def\mm{M} \def\ss{S}$Let \begin{align*} \f(x;\m,\s) &= \frac{1}{\sqrt{2\p}\s} e^{-(x-\m)^2/(2\s^2)} \end{align*} and \begin{align*} \f_m(x) &= \prod_{i=1}^m \f(x;\m_i,\s_i) \\ &= \frac{1}{(2\p)^{m/2}\prod_{i=1}^m \s_i} \exp\left(-\sum_{i=1}^m \frac{(x-\m_i)^2}{2\s_i^2}\right). \end{align*} By completing the square one finds \begin{align*} \f_m(x) &= A(m)\f(x;\mm,\ss), \end{align*} where \begin{align*} \frac{1}{\ss^2} &= \sum_{i=1}^m \frac{1}{\s_i^2} \\ \mm &= \sum_{i=1}^m \frac{\m_i}{\s_i^2} \\ A(m) &= \frac{\ss}{(2\p)^{(m-1)/2} \prod_{i=1}^m \s_i} \exp \left[ \frac12\left( \frac{\mm^2}{\ss^2} - \sum_{i=1}^m \frac{\m_i^2}{\s_i^2} \right) \right]. \end{align*} That is, such a product is itself proportional to a normal distribution. Thus, \begin{align*} \int_{-\infty}^\infty x^n \f_m(x) \,dx &= A(m) \int_{-\infty}^\infty x^n \f(x;\mm,\ss)\,dx \\ &= A(m) \int_{-\infty}^\infty ((x-\mm)+\mm)^n \f(x;\mm,\ss)\,dx \\ &= A(m) \int_{-\infty}^\infty \sum_{k=0}^n \binom{n}{k} (x-\mm)^k \mm^{n-k} \f(x;\mm,\ss)\,dx \\ &= A(m) \sum_{k=0}^n \binom{n}{k} \mm^{n-k} \int_{-\infty}^\infty (x-\mm)^k \f(x;\mm,\ss)\,dx \\ &= A(m)\sum_{k=0\atop k{\textrm{ even}}}^n \binom{n}{k} \mm^{n-k} \ss^k (k-1)!!, \end{align*} and so \begin{align*} \int_{-\infty}^\infty x^n \f_m(x) \,dx &= \frac{\ss \mm^n}{\prod_{i=1}^m \s_i} \frac{1}{(2\p)^{(m-1)/2}} \exp \left[ \frac 12\left( \frac{\mm^2}{\ss^2} - \sum_{i=1}^m \frac{\m_i^2}{\s_i^2} \right) \right] \\ & \quad \times \sum_{k=0}^{\lfloor n/2\rfloor} \binom{n}{2k} (2k-1)!! \left(\frac{\ss}{\mm}\right)^{2k} \end{align*}
For the original problem, \begin{align*} m &= 2 \\ (\m_1,\s_1) &= (0,\s) \\ (\m_2,\s_2) &= (\m,\r) \end{align*} so \begin{align*} \mm &= \frac{\m\s^2}{\s^2+\r^2} \\ \ss^2 &= \frac{\s^2\r^2}{\s^2+\r^2}. \end{align*} With a little work, we find a formula agreeing with that of @CarloBeenakker.
Yes, there are closed form expressions for $$I_n = \int_{-\infty}^{\infty} \frac{1}{\sqrt{2 \pi} \sigma} e^{-\frac{x^2}{2 \sigma^2}} \frac{1}{\sqrt{2 \pi} \rho} e^{-\frac{(x-\mu)^2}{2 \rho^2}} x^n\mathrm{d}x,\;\;n\in\mathbb{N}.$$
For $\mu=0$ this is simplest, one has $$I_n=(2\pi)^{-1/2}(n-1)\text{!!} \left(\rho^2 \sigma^2\right)^{n/2} \left(\rho^2+\sigma^2\right)^{-n/2-1/2},\;\;\text{for}\;\;n\;\;\text{even},$$ $I_n=0$ for $n$ odd.
For $\mu\neq 0$ the expressions are a bit more lengthy, $$I_n=Z_nP_n\left(\frac{\mu \sigma}{\rho \sqrt{\rho^2+\sigma^2}}\right),$$ $$Z_n=(2\pi)^{-1/2}e^{-\frac{\mu^2}{2 \left(\rho^2+\sigma^2\right)}} (\rho \sigma)^n \left(\rho^2+\sigma^2\right)^{-n/2-1/2}.$$ The function $P_n(x)$ is a polynomial of degree $n$, the first eight are given by $$\left\{x,x^2+1,x^3+3 x,x^4+6 x^2+3,x^5+10 x^3+15 x,x^6+15 x^4+45 x^2+15,x^7+21 x^5+105 x^3+105 x,x^8+28 x^6+210 x^4+420 x^2+105\right\}.$$ This is a polynomial described at https://oeis.org/A100861, given by $$P_n(x)=x^n+\sum _{k=1}^{n/2} \frac{n! x^{n-2 k}}{2^k k! (n-2 k)!}.$$
