Here, I have a recurrence relation:

F_{n+1}(ω) = (ω – α_{n})×F_{n}(ω) – β_{n-1}×F_{n-1}(ω),

with (obtaining these two coefficients below is my main goal)

α_{n} = [∫ω×F^2_{n}(ω)×h^2(ω)dω]/[∫F^2_{n}(ω)×h^2(ω)dω],

β_{n} = [∫F^2_{n}(ω)×h^2(ω)dω]/[∫F^2_{n-1}(ω)×h^2(ω)dω],

where the lower (upper) limit of integration is 0 (Λ).

The initial conditions satisfy F_{n<0}(ω)=0 and F_{0}(ω)=1. In fact, I wrote a Julia code to implement it, but due to the excessive function iterations, the computation becomes too slow for larger n. Therefore, I want to know if there is a more efficient method to obtain the coefficients α_{n} and β_{n}.

using QuadGK;


# initial parameters
Γ0 = 1e-2;
Λ = 20;


h_squared(ω) = (Γ0 / π^2) * (ω / (1 + ω)^2);


function recursion_coff(πk, pnπk)

    nπk = quadgk(ω -> πk(ω)^2 * h_squared(ω), 0, Λ)[1]

    αk = quadgk(ω -> ω * πk(ω)^2 * h_squared(ω), 0, Λ)[1] / nπk

    βk = nπk / pnπk
    return αk, βk, nπk
end;


N = 60;
em = [];
tm = [];

πkm1 = ω -> 0.0;
πk = ω -> 1.0;
pnπk = 1.0;


for k in 1:20
    αk, βk, pnπk = recursion_coff(πk, pnπk)

    πkp1 = let πk = πk, πkm1 = πkm1, αk = αk, βk = βk
        ω -> (ω - αk) * πk(ω) - βk * πkm1(ω)
    end

    πkm1 = πk
    πk = πkp1
    
    push!(em, αk)
    push!(tm, βk)

    println("Iteration: $k, αk = $αk, βk = $βk")
end

I want to know if it is possible to improve the efficiency of computing the recurrence relation of the function at the code level, not limited to using Python, Matlab, or Julia. Thanks everyone.