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

Question

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.