Simpson's method

#include <iostream>
#include <cmath>

using namespace std;

double f(double x) {
    return sqrt(1 + pow(x, 3));
}

double oddSum(double f(double x), int n, double h, double a) {
    double sum = 0;
    for (int i = 1; i < (n / 2); i++) sum += f(a + (2 * i - 1) * h);
    return sum;
}

double evenSum(double f(double x), int n, double h, double a) {
    double sum = 0;
    for (int i = 2; i < (n / 2); i++) sum += f(a + (2 * i * h));
    return sum;
}

double SimpsonsIntegral(double f(double x), double a, double b, double eps) {
    double intgrl = 0, prevIntgrl = 0, n = 2, h = (b - a) / n;

    do {
        prevIntgrl = intgrl;
        intgrl = (h * (f(a) + f(b) + 4 * evenSum(f, n, h, a) + 2 * oddSum(f, n, h, a))) / 3;

        n *= 2; h /=2;
    } while(fabs(prevIntgrl - intgrl) >= eps);

    return intgrl;
}

int main() {
    cout << "Simpson's method:" << endl;
    cout << "Integral = " << SimpsonsIntegral(f, 1, 3, 1e-4);

    return 0;
}