
% FILE doubleint0.m begins.
% Compute \int_{x_1}{x_2}  \int_{y_1}^{y_2} myf(x,y) dx dy
% This is based on the Riemann sum with m (n) steps in x (y) direction
function val=doubleint0(x1,x2,y1,y2,m,n)
a=x1; b=x2;
s=0.0;
for ii=0:(m-1)
      x=a+ii*(b-a)/m;
for jj=0:(n-1)
      dy=(y2-y1)/n;
      y=y1+jj*dy;
      s=s+myf(x,y)*dy*(b-a)/m;
end
end
val=s;
% fprintf('%12.5f\n', s)

function z= myf(x,y)
      z= x.*y;
% FILE doubleint0.m ends.
