#include <iostream>
#include <cmath>
#include <vector>

using namespace std;

// Function of given equation
double function(double x)
{
    double result;
        result = 0.2+25*x-200*pow(x,2)+675*pow(x,3)-900*pow(x,4)+400*pow(x,5);
    return result;
}

int main()

{
    double boundup;         // Upper bound
    double bounddown;       // Lower bound
    double n;               // Number of intervals
    double step;            // Step size
    double fsimp;           // Simpson's rule
    double ftrap;           // Trapezoidal rule

    cout << "Enter lower bound: " ;
    cin >> bounddown;

    cout << "Enter upper bound: " ;
    cin >> boundup;

    cout << "Number of intervals: " ;
    cin >> n;

    // Vector class for both x and fx
    vector<double> x(n+1);
    vector<double> fx(n+1);

    // Define first and last values in x and fx because they will not change
    step = (boundup-bounddown)/n;
    x[0] = bounddown;
    x[n] = boundup;

    fx[0] = function(x[0]);
    fx[n] = function(x[n]);

    fsimp = fx[0];

    for (int i=0; i<n; i++)
    {
        x[i+1] = x[i]+step;
        fx[i+1] = function(x[i+1]);

        if (i%2=0)
        {
            fsimp = 2*fx[i+1] + fsimp;
        }

        if (i%2=1)
        {
            fsimp = 4*fx[i+1] + fsimp;
        }

        fsimp = fx[n] + fsimp;
    }

    cout << "Bounds of integration: " << bounddown << ", " << boundup << endl;
    cout << "Number of intervals: " << n << 'n' << endl;
    cout << "Integral value of f(x): " << fsimp << endl;

return 0;
}