import numpy as np

def f(x):
    return np.exp(-((x-3)**2)/2)+2*np.exp(-((x-8)**2)/3)+0.5*np.exp(-((x-13)**2)/6)

def simpsonIntegral(f, low, high, N):
    step = (high - low) / (2*N)
    x0 = low + 3 * step
    x1 = low + 2 * step
    S0, Se = 0, 0
    for i in range(0, N - 1):
        S0 += f(x0)
        Se += f(x1)
        x0 += 2*step
        x1 += 2*step
    S0 += f(low + step)
    S0 *= (4 / 3) * step
    Se *= (2 / 3) * step
    sum = S0 + Se + ((step / 3) * (f(low) + f(high)))
    return sum

sum = simpsonIntegral(f,0,20,10)
print(f'The integral between {low} and {high} = {sum:.6f}')