from sympy import *
from sympy.parsing.sympy_parser import parse_expr

x, y, n = symbols('x y n')

def simpsons_rule(f, lower_bound, upper_bound, num_intervals):
    '''
    returns the approximation of an integral
    f = the function whos integral is to be approximated
    lower_bound = the lower bound of the integral
    upper_bound = the upper bound of the integral
    num_intervals = the number of intervals to calculate underneath the curve
    '''
    approximation = f.subs(x, lower_bound)
    delta_x = (upper_bound - lower_bound) / num_intervals
    this_x_value = lower_bound + delta_x
    #the second coefficient of the formula is 4
    coefficient = 4
    for _ in range(num_intervals - 1):
        approximation = approximation + coefficient * f.subs(x, this_x_value)
        this_x_value = this_x_value + delta_x
        #the coefficients of the formula then alternate between 2 and 4 except the last one
        if coefficient == 4:
            coefficient = 2
        else:
            coefficient = 4
    return (delta_x / 3 * (approximation + f.subs(x, upper_bound))).evalf()