% matlab script to evaluate two-dimensional integral
% once via Monte Carlo
% and once using integral2

close all
clear all
format long

% first Monte Carlo

f = @(x,y) 1./(1-x.*y);

Npts = input('Enter number of points per dimension: ');

sum = 0;
area = 1;

for i=1:Npts,
   for j=1:Npts,
      x = rand;
      y = rand;
      sum = sum + feval(f,x,y);
   end
end
MC = sum/(Npts^2)*area
exact = pi^2/6
errorMC = exact - MC

% now integral2

Q1 = integral2(f,0,1,0,1) % warnings!

errorQ1 = exact - Q1

Q2 = integral2(@(x,y) (1./(1-x.*y)),0,1-eps,0,1-eps)
% now we know to what extent we can trust the first answer

errorQ2 = exact - Q2