# R code to integrate h() using classical Monte Carlo method
# Example 2

ex2.fcn <- function(x){
h <- 4 - (x[1])^2 - 2*(x[2])^2
return(h)
}

my.n <- 10000

a1 <- 0; b1 <- 1
a2 <- 0; b2 <- 1

# If the range of integration over x1 is different from
# the range of integration over x2, then you would need to
# specify an [a1,b1] corresponding to x1 and 
# an [a2,b2] corresponding to x2.
# Here the range is [0,1] for both x1 and x2.

rand.unif <- cbind(runif(my.n, min=a1, max=b1), runif(my.n, min=a2,max=b2))
h.values <- apply(rand.unif, 1, ex2.fcn)

I.approx <- ((b1-a1)*(b2-a2)/my.n)*sum(h.values)

print(I.approx)