Why is my Simpson's modelling not giving the correct answer?

Question

I've been trying to do some Simpson's modelling to approximate integrals and although my code runs fine, it doesn't seem to give the correct answer. I'm new-ish to coding so I understand it's probably a dumb mistake but I can't figure out what's wrong. f(x) has already been defined and seems to work, it's just this integral approximation that's a bit dodgey.

def integrate_numeric(xmin, xmax, N):
    ''' 
    Numerical integral of f from xmin to xmax using Simpson's rule with 
        N panels.
    '''
    h = (xmax - xmin)/2*N
   
    # Odd sum
    oddtotal = 0
    for i in range(1, N, 2):
        oddtotal += f(xmin + (i*h))
    f_odd = oddtotal
    
    # Even sum
    eventotal = 0
    for i in range(2, N, 2):
        eventotal += f(xmin + (i*h))
    f_even = eventotal
    
    return (h/3)*(f(xmin)+(4*f_odd)+(2*f_even)+f(xmax))

Answer 1

• The ranges for odd numbers need to be from x=1 to x=N-1
• The ranges for even numbers need to be from x=2 to x=N-2
• h is needs to be (xmax - xmin)/N

Here is the Simpsons rule, with those corrections:

def f(x):
    return x**2


def integrate_numeric(xmin, xmax, N):
    '''
    Numerical integral of f from xmin to xmax using Simpson's rule with
        N panels.
    '''
    h = (xmax - xmin) / N

    odd_sum = sum(f(xmin + i*h) for i in range(1, N, 2))
    even_sum = sum(f(xmin + i*h) for i in range(2, N-1, 2))

    return h/3 * (f(xmin) + 4*odd_sum + 2*even_sum + f(xmax))


print(integrate_numeric(0, 10, 50)) -> 333.3333333333333

Alternatively, just use scipy's simpsons rule