[Math] Using Monte Carlo method to estimate the value of $\int_0^1\int_0^1e^{xy} \ dx \ dy$

Question

Monte Carlo, estimate the value of $$\int_0^1\int_0^1e^{xy} \ dx \ dy.$$

I am using matlab to solve this problem. My code is the following:

A = 1;  
N = 100000;  
s = 0;  
for i = 1:N;  
    x = rand; 
    y = rand;  
    z = rand;  
    if z <= exp(x*y);  
        s = s + 1;  
    end;  
end;  
I=A*s/N

The problem is that every time I run this script I get an answer of I = $1$. In fact, whatever I initialize the area (in this case volume) to be is what gets returned as the answer. As a visual aid, I attempted to graph $e^{xy}$.

Now, this code is an adaptation of a code I wrote for a previous problem I was working on. That problem was to use the Monte Carlo method to integrate $\int_{0}^{1}x^2dx$ in matlab. My code for that problem was the following:

A=1; 
N=10000; 
s=0; 
for i=1:N 
    x=rand; 
    y=rand; 
        if y<= x^2; 
            s=s+1; 
        end; 
end; 
I=A*s/N

This code was successful, and, as I increased the value for N, it became increasingly more accurate, as expected. However, when I make the jump from the single integral to the double integral, something in my code is not translating properly.

Answer 1

Here is a code that should do the trick:

s=0;
for i=1:N  
x=rand(0,1); 
y=rand(0,1);
s=s+exp(xy); 
end 
I = s/N;