I have tried to solve the following integral in closed form: (I encountered this on approximating a CRLB)
\begin{equation} \int_{-\infty}^{\infty} x^2 \phi(x) \phi(a+bx) dx \end{equation}
where $\phi(x)=\frac{1}{\sqrt{2 \pi}} e^{-\frac{x^2}{2}}$ is the pdf of a standard Gaussian Random Variable. Using the 'integration by parts' and Owen's 'a table of normal integrals' I found following solution, but integrating the function Numerically (In Matlab) gives me different values. Can any body help me? There must be problem with my closed-form solution
\begin{equation} \int_{-\infty}^{-\infty} x^2 \phi(x) \phi(a+bx) dx= \frac{1}{(1+b^2)^{3/2}} ~ \phi\left(\frac{a}{\sqrt{1+b^2}}\right) \end{equation}