Trapezoidal Rule Code

I am including the code I wrote for the trapezoid rule integration assignment in this post. In addition to the original assignment, I incorporated an example of an array, per request. My original intent for usage of arrays fell through, and so I am attempting a weak make-up array display. For the sake of the reader, I will include comments besides a few of the lines.


#include <stdio.h>
#include <cmath>


double trapezoid (int n, float a, float b);             //Declaring both functions which will be utilized later
double f (float x);

int main() {

        float variable[2] = { 2, 6 };           //Creation of an array name "variable", with type float
        float a = variable[0];                  //Both a and b are assigned values of elements from the array
        float b = variable[1];                  //Array elements start at "0", so the array of two goes "0,1"
        double val;
        double answer = 32.0/3.0;
        double minErr = 1000;                   //These lines of code will be used to get an optimal n-value
        double err;
        int optimalN;
        double optimalVal;

        for( int n = 0; n < 5000; n++ ) {       //Creation of a for loop to work the integration

        val = trapezoid (n, a, b);

        err = std::abs(answer-val);

        if( err < minErr ) {                    //If statement that is used to capture and store the optimal values
            minErr = err;
            optimalN = n;
            optimalVal = val;
          }
        }

        printf ("The integral value = %f\n", val);
        printf( "minimum err = %f\n", minErr );
        printf( "Optimal n = %d\n", optimalN );
        printf( "Optimal val = %f\n", optimalVal );


        return 0;

}
                                                //It is considered good programming to define functions at the end- "prototyping"
double f (float x) {                            //Here, the first function is defined, with full parameters and statements

        return 4.0-(x-4.0)*(x-4.0);
        }

double trapezoid (int n, float a, float b) {    //The second function does all the integration

        double sum = 0;
        float step = (b-a)/n;
        float xi = a;

        for (int i = 0; i < n; i++) {       //This for loop integrates for each trapezoid of n, up to the max 5000
         sum += step*((f(xi+step)+f(xi))/2);
         xi += step;
         }
        return sum;

        }