double simpson(double(*f)(double), double a, double b, int n)
/*
    Integration by integration by Simpson rule of f(x) on [a,b]
    method:     Simpson rule
    written by: Alex Godunov
-------------------------------------------------
    f - Function to integrate (supplied by a user)
    a - Lower limit of integration
    b - Upper limit of integration
    s - Result of integration (out)
    n - number of intervals
=================================================
*/
{
    double s, dx, x;
// if n is odd we add +1 interval to make it even
    if((n/2)*2 != n) {n=n+1;}
    s = 0.0;
    dx = (b-a)/static_cast<float>(n);
    for ( int i=1; i<=n-1; i=i+2)
    {
        x = a+static_cast<float>(i)*dx;
        s = s + 4.0*f(x) + 2.0*f(x+dx);
    }
    s = (s + f(a)+f(b))*dx/3.0;
    return s;
}