def simpson(f, a, b, n):
    """f=function, a=initial value, b=end value, n=number of double intervals of size 2h"""

    n *= 2
    h = (b - a) / n;
    S = f(a)

    for i in range(1, n, 2):
        x = a + h * i
        S += 4 * f(x)

    for i in range(2, n-1, 2):
        x = a + h * i
        S += 2 * f(x)

    S += f(b)
    F = h * S / 3

    return F