Integral approximation using gauss-hermit quadrature method: mvQuad package

You should use the product of PDFs for both X and Y in function myFun2d

myFun2d <- function(x) {
  x[, 1] * exp(-0.5 * (x[, 1]-0.5)^2) / sqrt(2 * pi) * x[, 2] * exp(-0.5 * (x[, 2]-0.5)^2) / sqrt(2 * pi)
}

and you will see

> library(mvQuad)

> # create grid
> nw <- createNIGrid(dim = 2, type = "GHe", level = c(10, 10))

> m <- c(0.5, 0.5)

> c <- matrix(c(1, 0, 0, 1), nrow = 2, byrow = F)

> rescale(nw, m = m, C = c, dec.type = 0)

> # define the integrand
> myFun2d <- function(x) {
+   x[, 1] * exp(-0.5 * (x[, 1] - 0.5)^2) / sqrt(2 * pi) * x[, 2] * exp(-0.5 * (x[, 2] - 0.5)^2) / .... [TRUNCATED]

> # compute the approximated value of the integral
> (A <- quadrature(myFun2d, grid = nw))
[1] 0.25