I have to do a double integral of the function sin(x) + cos(y) dxdy in IML , using a Monte Carlo simulation, with vectorized code. Truly, I do not know how to do it.
There are two common ways to use Monte Carlo to estimate integrals, so you'll have to look at your notes/textbook to decide which you should use. In one (the average-value method), you generate points uniformly on a 2-D domain and evaluate the function at each random location. In the other method (the volume method), you generate 3-D points and estimate the probability that a random point is under the graph of the function.
Here are some articles that use the average-value method. Study and understand the 1-D example before you attempt 2-D.