Оцените двойные интегралы, используя составное правило Симпсонов

Ранее я создал этот код для оценки интеграла с использованием составного правила Симпсона, и теперь мне нужно изменить его для оценки двойных интегралов. У кого-нибудь есть идеи, как это сделать? Буду благодарен за любую помощь ..

#include <iostream>
#include <cmath>

using namespace std;

float f(float x); //returns the function we are integrating
void simpsons_rule(float a, float b, int n); // Carries out Simpsons composite ruleint main()
{
cout << "This program finds an approximation of the integral of a function under two limits using Simpsons composite rule." << endl;
// User is introduced to program and given explanation of what it does
float a,b;
int n;
cout << "Enter the lower limit " << endl;
cin >> a;
cout << "Enter the upper limit " << endl;
cin >> b; // User has freedom to specify the upper and lower limits under which the function will be integrated
cout << "Enter the number of intervals (must be an even number) " << endl;

do
{
cin >> n; // User can also choose the number of intervals

if(n%2 == 0 && n!=0) // If an even number of intervals was entered, the program continues and simpsons rule is carried out
{
simpsons_rule(a,b,n);
}
else
{
cout << "Invalid number of intervals, please re-enter a value" << endl; //If an odd number of intervals was entered, the user is prompted to enter a valid number
}
}while(cin); // do-while loop used to ensure even number of intervals was entered

return 0;
}

float f(float x)
{
return(pow(x,4) - 3*pow(x,3) + pow(x,2) + x + 1); // function to be integrated
}

void simpsons_rule(float a, float b, int n)
{
float s2 = 0;
float s3 = 0;
int i;
float h = (b-a)/n; // step length is calculated

for(i=1; i <= (n/2)-1; i++)
{

if((i%2) == 0) // checks if i is even
{
s2 = s2 + f(a+(i*h)); // Finds the sum of all the values of f(x) where x is even

}
}

for(i=1; i<=(n/2); i++)
{
if((i%2) == 1) // checks if i is odd
{
s3 = s3 + f(a+2*(i*h)); // Finds the sum of all the values of f(x) where x is odd
}
}

float s = (h/3)*(f(a)+ f(b) + (2*s2) + (4*s3)); // This variable contains the final approximaton of the integral

cout << "The value of the integral under the specified limits is: " << s << endl;