function DoubleIntegralRectangle(x_min, x_max, N_x, y_min, y_max, N_y,...
				 FunctionString)

% Computes the double integral of z=f(x,y) over the domain
% x_min <= x <= x_max, y_min <= y <= y_max, with N_x grid points
% in the x-direction and N_y grid points in the y-direction.
% The integral is computed using both the UpperRightCorner rule
% and the Midpoint rule
% The function is specified by the string FunctionString.
% For example if F(x,y)=x^2 + y^2 the value of FunctionString would be 
% 'X.^2 + 3*Y.^2'
%
% Uncomment the lines "Function(X,Y)" to see the values of f being summed
% Uncomment the call to GraphZvsXY to plot the surface z=f(x,y) 
%
% Typical call:
%
% DoubleIntegralRectangle(0,1,10,0,1,10,'3*x.^2+4*y.^2')

     Function = inline(FunctionString);

     % UpperRightCorner rule:

     delta_x = (x_max-x_min)/N_x;

     x = [x_min+delta_x:delta_x:x_max];

     delta_y = (y_max-y_min)/N_y;

     y = [y_min+delta_y:delta_y:y_max];

     [X,Y] = meshgrid(x,y);

     %Function(X,Y)

     UpperRightCornerDoubleSum = sum(sum(Function(X,Y)))*delta_x*delta_y

     % Midpoint rule:

     delta_x = (x_max-x_min)/N_x;

     x = [x_min+0.5*delta_x:delta_x:x_max-0.5*delta_x];

     delta_y = (y_max-y_min)/N_y;

     y = [y_min+0.5*delta_y:delta_y:y_max-0.5*delta_y];

     [X,Y] = meshgrid(x,y);

     %Function(X,Y)

     MidpointDoubleSum = sum(sum(Function(X,Y)))*delta_x*delta_y

     %GraphZvsXY(0,1,10,0,1,10,'3*x.^2+4*y.^2')