Hi,
You can refer to the code below.
function [Q] = gaussQuad2(f, xs, ys, n)
% Gauss quadrature for double integral over a rectangular domain
% n - number of integration points in each direction
% Define Gauss quadrature weights and nodes
[xs1, ws1] = gauss(n, xs(1), xs(3));
[xs2, ws2] = gauss(n, ys(1), ys(3));
% Compute integral approximation
Q = 0;
for i = 1:n
for j = 1:n
Q = Q + ws1(i)*ws2(j)*f(xs1(i),xs2(j));
end
end
f = @(x,y) x*y;
xs = [0, 2, 4, 0];
ys = [0, 0, 4, 4];
n = 4;
Q = gaussQuad2(f, xs, ys, n);
disp(Q);
