# Simpson's Rule to approximate integral( f(x), x, a, b ) over N (even) subintervals
# The user input is between the hashtags. You can run the code by cutting and 
# pasting into a SageMath Cell.
#
def f(x):
    return x**2 + x
a = 0.0
b = 3.0
N = 30  # Must be even
#
h = ( b - a ) / (1.0*N)
h2 = 2 * h
s0 = f(a) + f(b)
x = a - h
s1 = 0.0
for i in range(0,N,2):
    x = x + h2
    s1 = s1 + f(x)
s1 = 4 * s1
x = float(a)
s2 = 0.0
for i in range(1,N-1,2):
    x = x + h2
    s2 = s2 + f(x)
s2 = 2 * s2
s = h * ( s0 + s1 + s2 ) / 3.0
print( "Simpson's rule approximation:", s )