Implementing Simpson's rule with two variables

For one of my lab project I need to integrate a given function (in this example, f1, using Simpson's ruls. When only one variable x exists, I was able to find a code that calculate the integration:

f1[x_] := 3 + Cos[x + 3]
LEFT1[f_, a_, b_, n_] :=
Module[{h, i},
h = (b - a)/n;
N[h Sum[f1[a + i h], {i, 1, n - 1, 2}]]]
RIGHT1[f_, a_, b_, n_] :=
Module[{h, i},
h = (b - a)/n;
N[h Sum[f[a + i h], {i, 2, n - 2, 2}]]]
TRAP1[f_, a_, b_, n_] := 4 LEFT1[f, a, b, n] + 2 RIGHT1[f, a, b, n]
MID1[f_, a_, b_, n_] :=
Module[{h, i},
h = (b - a)/n;
N[h (f[a] + f[b])]]
SIMP1[f_, a_, b_, n_] := (MID1[f, a, b, n] + TRAP1[f, a, b, n])/3.0

SIMP1[f1, -1, 10, 10]

32.5062

Instead of using f1[x_] := 3 + Cos[x + 3], I would like to use f1[x_] = 3 + Cos[x + 3 b] and calculate the integration for value of b ranging from 1 to 50. I am not sure how I should start tackling this problem.

How about this? It implements the tensor product rule associated to the quadrature weights of Simpson 1D-rule.

Simpson2D[f_, {xa_, xb_, xn_}, {ya_, yb_, yn_}] :=
Module[{Δx, xω, Δy, yω},
Δx = (xb - xa)/(2. xn);
xω = ConstantArray[2./3. Δx, 2 xn + 1];
xω[[2 ;; ;; 2]] = 4./3. Δx;
xω[] = 1./3. Δx;
xω[[-1]] = 1./3. Δx;

Δy = (yb - ya)/(2. yn);
yω = ConstantArray[2./3. Δy, 2 yn + 1];
yω[[2 ;; ;; 2]] = 4./3. Δy;
yω[] = 1./3. Δy;
yω[[-1]] = 1./3. Δy;
xω.Outer[f, Subdivide[N[xa], N[xb], 2 xn], Subdivide[N[ya], N[yb], 2 yn]].yω
]

Test:

f = {x, y} \[Function] 3. + Cos[x + y];
a = Simpson2D[f, {-1, 10, 100}, {1, 50, 100}];
b = Integrate[f[x, y], {x, -1, 10}, {y, 1, 50}];
Abs[1 - a/b]

3.21721*10^-9

• Thanks Henrik. I may have misstated my problem. I would need to print the integration over x for a specific value of b. For example, let assume that 1<b<10 with incrementation of 0.2. I would need to print the integration over x for b = 1,1.2,1.4,.1.6....10. – Jean-Hubert Olivier Aug 8 at 14:05