Consider the function $f(x_1,x_2,x_3) = \frac{1}{\sqrt{2\pi^3}}*\Bigl(e^{\frac{1}{2}(x_1^2+x_2^2+x_3^2)}\Bigr)*\Bigl(1+x_1x_2x_3*e^{\frac{1}{2}(x_1^2+x_2^2+x_3^2)}\Bigr)$. I am trying to find its marginal density with respect to the vector $\textbf{X} = (x_1,x_2)^{T}$. The worked example states that the result should be equal to $f(x_1,x_2)=\frac{1}{2\pi}*e^{-\frac{1}{2}(x_1^2+x_2^2)}$ because the function $g(x_3)=x_1x_2x_3*e^{\frac{1}{2}(x_1^2+x_2^2+x_3^2)}$ is odd, so when you integrate it over the interval $(-\infty,\infty)$ you get $\int_{-\infty}^\infty g(x_3) dx_3 = 0$. However, when I attempt it myself I fail in recreating the solution given in the example.
My progress so far
I have tried integration by substitution but that really didn't get me very far. Call $u(x_3) = \Bigl(e^{\frac{1}{2}(x_1^2+x_2^2+x_3^2)}\Bigr)$, $v'(x_3) = \Bigl(1+x_1x_2x_3*e^{\frac{1}{2}(x_1^2+x_2^2+x_3^2)}\Bigr)$. Then it follows that $\int_{-\infty}^{\infty} f(x_3) dx_3 =\frac{1}{\sqrt{2\pi^3}}*\int_{-\infty}^{\infty} u(x_3)*v'(x_3) dx_3$.
Integration by substitution gives $\frac{1}{\sqrt{2\pi^3}}*\Bigl(e^{\frac{1}{2}(x_1^2+x_2^2+x_3^2)}\Bigr)*\int_{-\infty}^{\infty}\Bigl(1+x_1x_2x_3*e^{\frac{1}{2}(x_1^2+x_2^2+x_3^2)}\Bigr) dx_3 - \frac{1}{\sqrt{2\pi^3}}*\left(\int_{-\infty}^{\infty} {-x_3*\Bigl(e^{\frac{1}{2}(x_1^2+x_2^2+x_3^2)}\Bigr)}*\left(\int_{-\infty}^{\infty} 1+x_1x_2x_3*e^{\frac{1}{2}(x_1^2+x_2^2+x_3^2)}dx_3\right)dx_3\right)$.
From here on, however, I tried many different algebraic manipulations which all led me nowhere. So my question becomes - what do I do to get to the result given in the introduction?
Consider that $I_1(x_1,x_2):=\frac{1}{\pi^{1/2}}\int_{\mathbb{R}}x_1x_2x_3e^{-\|x\|^2}dx_3=x_1x_2(\int_\mathbb{R}x_3\frac{e^{-x_3^2}}{\pi^{1/2}}dx_3)e^{-(x_1^2+x_2^2)}=0$ because the integral in the parenthesis is the mean of a normal random variable $\sim \mathcal{N}(0,1/2)$. For a similar reason: $I_2(x_1,x_2):={(2\pi)^{-1/2}}\int_{\mathbb{R}}e^{-\|x\|^2/2}dx_3={(\int_\mathbb{R}{(2\pi)^{-1/2}}{e^{-x_3^2/2}}dx_3)}e^{-(x_1^2+x_2^2)/2}=e^{-(x_1^2+x_2^2)/2}$ because the integral in the parenthesis is the integral of the density of a normal random variable $\sim \mathcal{N}(0,1)$. So we have $$\begin{aligned}\int_\mathbb{R}f(x_1,x_2,x_3)dx_3&=\frac{1}{2^{3/2}\pi^{3/2}}((2\pi)^{1/2}I_2(x_1,x_2)+\pi^{1/2}I_1(x_1,x_2))\\ &=\frac{1}{2^{3/2}\pi^{3/2}}\bigg((2\pi)^{1/2}e^{-(x_1^2+x_2^2)/2}\bigg)\\ &=\frac{1}{2\pi}e^{-(x_1^2+x_2^2)/2} \end{aligned}$$