Интеграция с переменными пределами

Мне нужно оценить двойной интеграл, где внутренняя верхняя граница является переменной:
integral2 between -5 and 5 ( integral1 between 0 and y f(x)dx )dy,

Я застрял в расчете внешнего цикла, который зависит от внутреннего цикла. Мой код работает очень долго, но возвращает ноль.

Как я могу вычислить интеграл с переменными пределами?

Сначала я создал функцию doubleIntegrate, В первую очередь функция удерживает массивы с коэффициентами для правила трапеции.

double NumericIntegrationDouble::doubleIntegrate(double (*doubleFunc
(const double &x), double dy, const double &innerLowBound, const double
&outerLowBound)
{

double innerValue = 0.0;
double outerValue = 0.0;

// arrays which store function values for the inner (X) and the outer (Y) integration loop

// vector filled with coefficients for the inner poop (trapezoidal rule)
std::vector<double> vecCoeffsX(numberOfIntervalsDouble+1, 2);
vecCoeffsX[0] = 1;                         // fist coeff = 1
vecCoeffsX[vecCoeffsX.size()-1] = 1;       // last coeff = 1
std::vector<double> funcValuesX(numberOfIntervalsDouble+1);// vector filled with coefficients for the inner poop (trapezoidal rule)
std::vector<double> vecCoeffsY(numberOfIntervalsDouble+1, 2);
vecCoeffsY[0] = 1;                        // same as above
vecCoeffsY[vecCoeffsY.size()-1] = 1;      // same as above
std::vector<double> funcValuesY(numberOfIntervalsDouble+1)// Then i created a loop in a loop where dy and dy stands for step size of integration. The variables xi and yi stand for the current x and y value.

// outer integration loop dy
for(int i=0; i<=numberOfIntervalsDouble; i++)
{
double yi = outerLowBound + dy*i;
funcValuesY[i] = (*doubleFunc)(yi);

// inner integration loop dx
for(int j=0; j<=numberOfIntervalsDouble; j++)
{
double dx = abs(yi - innerLowBound) / (double)numberOfIntervalsDouble;
double xi = innerLowBound + j*dx;
funcValuesX[j] = (*doubleFunc)(xi);
double multValueX = std::inner_product(vecCoeffsX.begin(), vecCoeffsX.end(), funcValuesX.begin(), 0.0);
double innerValue = 0.5 * dx * multValueX;
suminnerValue = suminnerValue + innerValue;
}
//auto multValueY = std::inner_product(vecCoeffsY.begin(), vecCoeffsY.end(), funcValuesY.begin(), 0.0);

outerValue = 0.5 * dy * suminnerValue;
}
return outerValue;
}