Comecei a trabalhar neste livro (Exercício de Física Computacional 5.4) e seus exercícios e fiquei preso à seguinte pergunta:

Escreva uma função Python J (m, x) que calcule o valor de Jm (x) usando a regra de Simpson com N = 1000 pontos. Use sua função em um programa para fazer um gráfico, em um único gráfico, das funções de Bessel J0, J1 e J2 como uma função de x de x = 0 ex = 20.

Eu criei o código a seguir para avaliar a primeira parte da pergunta, mas não tenho certeza se isso está correto:

def f(x, t):
    return 1 / pi * (math.cos(x - t * math.sin(x)))

def float_range(a, b, c):
    while a < b:
        yield a
        a += c  
    N = 1000
    a = 0.0
    b = 20.0
    h = (b - a) / N
    c = 0.0
    d = pi
    h2 = (d - c) / N
    s = 0.5 * f(a, 1) + 0.5 * f(b, 1)
    s / 3
    S1 = 0
    S2 = 0
    for k in range(1, N):
        for j in range(0, N):
            if k%2 == 0:
                S1 += 2 / 3 * f(a + k * h, c + k * h2)
            else:
                S2 += 4 / 3 * f(a + k * h, c + k * h2)
    s += S1 + S2
    print(h * s)

Alguém poderia me ajudar a resolver esta questão, nunca usei funções de Bessel antes?