library(emulator)
Fun <- function(a,b) sqrt((1/(2*pi)))*exp(-0.5*(a-b)^2)
n <- 5000
x <- rnorm(n)
y <- rnorm(n)
u <- rnorm(n)
I <- quad.form(outer(x,x,Fun)*outer(y,y,Fun),u)
}