Don't know what is wrong with this double integral code

I'm trying to write a script that can calculate the double integral of a function. What I want to do is this:

• move along the y-axis.
• for each value of y, I calculate the area of the cross-section of the 3D graph, which is the integral of the function for x at that specific value of y.
• I multiply this area by a small increment in y to get an elemental volume.
• I add up all these volumes up to get the final volume under the graph, which is equal to the double integral.

I don't know why my code isn't giving me the write answer for sin(x)*sin(y) between 0 and pi for x limits and 0 and pi for y limits.

I attach my code as a .m file and also below:

f = @(a, b) sin(a)*sin(b);

x0 = 0;

xn = pi; % Integral limits

n = 1000*pi; % No. intervals for x

y0 = 0;

ym = pi; % Integral limits

m = 1000*pi; % No. intervals for y

hx = (xn - x0)/n; % x step size

hy = (ym - y0)/m; % y step size

I = 0;

for y = y0:hy:ym

A = 0;

for x = x0:hx:xn

A = A + hx*f(x,y);

A = A - 0.5*hx*(f(x0,y) + f(xn,y));

I = I + A*hy;

end

end

Answer 1

[Bjorn Gustavsson]

Use integral2 instead. That should solve this type of problem.

HTH