Несобственный интеграл

import numpy as np
 
def fun_x(x,xo):
    return np.exp(-(x*x*(x-xo)*(x-xo)/4))
 
eps = 0.00001 # требуемая точность
xo = -3
dx = 1
while xo < 3 + eps:
    a = -2 # первое приближение. Значение функции в этой точке < eps
    b = 2
    
    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 = fun_x(x,xo)
 
        #Ih1_s1 (Ih2_s1)) - Вспомогательная переменная, элемент формулы Симпсона
#Ih1_s2 (Ih2_s2)- Вспомогательная переменная, элемент формулы Симпсона
        Ih1_s1 = fx[1::2].sum()
        Ih1_s2 = fx[2:-1:2].sum()
        # Считаем интеграл с шагом h
        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()
        # Считаем интеграл с шагом 2h
        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:.3f}\nс точностью: {Ih_check:.3e}\n')
            break
        # если точность больше eps, удваиваем пределы
        a <<= 1
        b <<= 1
    xo += dx

import numpy as np
 
def fun_x(x,xo):
    return np.exp(-(x*x/4+(x-xo)*(x-xo)))
 
eps = 0.00001 # требуемая точность
xo = -3
dx = 1
while xo < 3 + eps:
    a = -4 # первое приближение. Значение функции в этой точке < eps
    b = 4
 
    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 = fun_x(x,xo)
 
        #Ih1_s1 (Ih2_s1)) - Вспомогательная переменная, элемент формулы Симпсона
#Ih1_s2 (Ih2_s2)- Вспомогательная переменная, элемент формулы Симпсtона
        Ih1_s1 = fx[1::2].sum()
        Ih1_s2 = fx[2:-1:2].sum()
        # Считаем интеграл с шагом h
        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()
        # Считаем интеграл с шагом 2h
        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}')
            print(f'control при X0={xo} -> {4*np.pi/5/np.exp(2*xo*xo/5):.3f}\n')
            break
        # если точность больше eps, удваиваем пределы
        a <<= 1
        b <<= 1
    xo += dx