// This function returns the result of using Simpson's rule over n subintervals
// to approximate the integral of f(x) from x=a to x=b.  
// n MUST be even---This is not checked is this function!  Check in main()!
//
// Notice that this function is being passed a function as an argument.
// double (*f)( double x ) is a pointer to a function that takes a 
// type-double argument and returns a double.  The () around *f are
// required.

double simpson( double (*f)( double x ), double a, double b, int n = 100 )
{
    int i;
    
    double h = ( b - a ) / n;
    double sum1 = 0.0, sum2 = 0.0, sum4 = 0.0;
    
    sum1 = (*f)(a) + (*f)(b);  // You could simply use f(a) and f(b)
    
    for ( i = 1; i <= n-1; i += 2 ) {
        sum4 = sum4 + (*f)( a + i*h );  // You could simply use f( a + i*h )
    }
    for ( i = 2; i <= n-2; i += 2 ) {
        sum2 = sum2 + (*f)( a + i*h );  // You could use f( a + i*h )
    }
    
    double approx = h / 3 * ( sum1 + 4.0 * sum4 + 2.0 * sum2 );
    
    return( approx );
}