# 2-dimensional Gauss-Hermite quadrature in R

A similar question was asked here and the given answer is perfect for a unidimensional integration.

I need to make bidimensional integration in R with a Gauss-Hermite quadrature: $$\int_{R^2} h(p1,p2) \phi(p1,p2) dp1 dp2$$ with $$p1$$ and $$p2$$ two parameters that follow a multivariate normal distribution of density $$\phi(.)$$ such as: $$p1,p2 \sim \mathcal{N_2} (\mu, \Sigma)$$ with $$\mu = \left(\matrix{0\\ 1} \right)$$ and $$\Sigma = \left[\matrix{0.6&0.5\\0.5&0.25} \right]$$.

I tried some packages in R but the results I got were very different from one another. I also tried the solution here but I'm not sure of the results I got.

I am looking for a simple solution, preferably based on gauss.quad.

Not being an expert, I have no idea of when to divide by $$\pi$$, if I should make a rescaling in my case and so on.

Thanks in advance for any help !