Simpson's method for approximating definite integral

simpson_method.py

# Numerical Analysis
# Simpson's Method
# Author: A.F. Agarap

# Function = (h / 3) (f(x[0]) + 4f(x[1]) + 2f(x[2]) + 4f(x[3]) + 2f(x[4]) + ... + 4f(x[n-1]) + f(x[n]))
# f(x) = cos(x**2) * exp(x**3)

import definite_integral as di
import os

def main():
	formula = input("Enter f(x): ")
	a = int(input("Enter value for a: "))
	b = int(input("Enter value for b: "))
	n = int(input("Enter value for n: "))
	while (n % 2 != 0):
		n = int(input("Enter value for n: "))
	x = di.generate_set(a, b, n)
	result = 0
	for i in range(0, (n + 1)):
		if (i == 0) or (i == n):
			result += di.function(formula, x[i])
		elif (i % 2 == 0):
			result += (2 * di.function(formula, x[i]))
		elif (i % 2 != 0 and (i != 0 or i != n)):
			result += (4 * di.function(formula, x[i]))
	result *= ((abs(a - b) / n) / 3)
	print(result)
	# os.system('pause') # For Windows OS

main()

The source program for definite_integral:

from math import *

def generate_set(a, b, n):
x = []
h = abs(a - b) / n
for i in range(n + 1):
x.append(a + (i \* h))
return x

def function(formula, x):
# return math.cos(x**2) \* math.exp(x**3)
return eval(formula)