Решение интеграла через метод Симпсона

import numpy as np
 
def fun1(x,x0):
    return np.exp(2*(np.sqrt(2)*((np.exp(-(x*x)))/(((x0*x)*(x0*x)/9)+1/9))))
 
def fun2(x,x0):
    return np.exp((np.sqrt(2)*np.sqrt(np.pi))/2)
    
def fun3(x,x0):
    return np.exp((81*np.pi)/2)
 
def integr(f, eps=1e-5, xo=0):
    a = -8
    b = 8
 
    while True:
        n = int(np.ceil((b-a)/(eps**(0.25))/4)*4) 
        n = max(n, 40)
        #print('n =', n)
        h = (b-a)/n 
        x = np.linspace(a, b, n+1) 
        fx = f(x,xo)
 
        Ih1_s1 = fx[1::2].sum()
        Ih1_s2 = fx[2:-1:2].sum()
        
        Ih1 = h/3*(fx[0]+fx[40]+4*(Ih1_s1)+2*(Ih1_s2))  
 
        Ih2_s1 = fx[2:-1:4].sum()
        Ih2_s2 = fx[4:-4:4].sum()
        
        Ih2 = 2*h/3*(fx[0]+fx[40]+4*(Ih2_s1)+2*(Ih2_s2))
        Ih_check = abs(Ih1-Ih2)/15 
 
        if Ih_check < eps:
           # print(f'Значение интеграла при X0={xo} -> {Ih2**2:.3f}\nс точностью: {Ih_check:.3e}')
            return Ih2
            
        a <<= 1
        b <<= 1
    xo += dx
 
x1 = -6
x2 = 6
dx = 1
Int_2 = integr(fun2, 1e-5)
Int_3 = integr(fun3, 1e-5)
 
while x1<=x2:
    Int_1 = integr(fun1, 1e-5, x1)
    print(f'При xo={x1} отношение интегралов равно -> {Int_1**2/(Int_2*Int_3):.5f}')
    x1 += dx