Numerical solution of integro-differential equation where the integral term contains a double integral

I have the following code:

sys = {
D[f[ϵ, x], {x, 2}] + Iee[ϵ, x] + Iph[ϵ, x] == 0,
f[ϵ, -L/2] == fF[ϵ - eV/2],
f[ϵ, L/2] == fF[ϵ + eV/2]
}
where:
Iee[ϵ_, x_] := a Integrate[
Integrate[
(f[ϵ, x] f[ϵ1 - ω, x] (1 - f[ϵ - ω, x]) (1 - f[ϵ1, x]) -
f[ϵ - ω, x] f[ϵ1, x] (1 - f[ϵ, x]) (1 - f[ϵ1 - ω, x])),
{ϵ1, -∞, ∞},
{ω, -∞, ∞}
]
]

Iph[ϵ_, x_] := b Integrate[
ω^2 ((1 - f[ϵ, x]) f[ϵ + ω, x] (1 + n[ω]) +
(1 - f[ϵ, x]) f[ϵ - ω] n[ω] -
f[ϵ, x] (1 - f[ϵ - ω]) (1 + n[ω]) -
f[ϵ, x] (1 - f[ϵ + ω]) n[ω]),
{ω, 0, ∞}
]
n[ω_] = 1/(Exp[ω/T] - 1);
fF[ϵ_] = 1/(1 + Exp[ϵ/T]);


where a,b,T,L,eV are constants.

I need to solve this equation numerically. As I understand, simple functions as well as NDSolve don't work here, and I need to use The Numerical Method of Lines (like here Integrating over $x$ in numerically solving a partial integrodifferential equation)

But anyway I still dont understand what to do with limits of integration and how to connect it with my boundary conditions.

• For a numerical solution it would be helpful to find finite integration limits. I don't know the parameter values Do you think it's possible to find an appropriate transformation w,ϵ? Dec 13, 2023 at 20:03
• As the first step towards a numerical solution, it would be necessary to determine the parameters a,b,T, L, eV. Dec 14, 2023 at 3:35
• @UlrichNeumann Im not sure, but i gonna try to use something like 2 last functions for such transform. It might be naturally in this contecst. Dec 14, 2023 at 22:37
• @AlexTrounev I thought that its possible to solve for all real. For the first solution it should be ok use any non zero constants. Dec 14, 2023 at 22:42
• @PabloEscobar In Iph definition there are f[ϵ, x] and f[ϵ - ω]. Is it a typo? Dec 15, 2023 at 14:13