double trapezoid(double(*f)(double), double a, double b, int n)
/*
    Integration by trapezoid rule of f(x) on [a,b]
    method:     trapezoid rule
    written by: Alex Godunov
-------------------------------------------------
    f - Function to integrate (supplied by a user)
    a - Lower limit of integration
    b - Upper limit of integration
    r - Result of integration (out)
    n - number of intervals
=================================================
*/
{
    double r, dx, x;
    r = 0.0;
    dx = (b-a)/static_cast<float>(n);
    for (int i = 1; i <= n-1; i=i+1)
    {
        x = a + static_cast<float>(i)*dx;
        r = r + f(x);
    }
    r = (r + (f(a)+f(b)/2.0))*dx;
    return r;
}