function J = simpsondouble(f, c, d, a, b, m, n)
% Integrate f(x,y) over c(x) <= y <= d(x) using Simpson's rule.
% The function handles must accept array inputs.

x = linspace(a, b, n + 1);
dy = (d(x) - c(x)) / m;
y = c(x) + (0:m)' .* dy;

wx = [1, 2 + 2 * mod(1:n - 1, 2), 1];
wy = [1; 2 + 2 * mod((1:m - 1)', 2); 1];

inner_integral = dy / 3 .* sum(wy .* f(x, y), 1);
J = (b - a) / (3 * n) * sum(wx .* inner_integral);
