Incorrect integration results depending on integration bounds, using quadgk directly works

[cpf-work]

This comes from obtaining electron densities in a density of states given a chemical potential eta. Here's a MWE which shows the issue

using Integrals
using QuadGK
using Plots

function integrand(E::Real, eta::Real)
    return exp(-E^2/2) / (sqrt(2π) * (exp(E-eta) + 1))
end

function using_quadgk(eta::Real, lim::Real)
    iq(x) = integrand(x, eta)
    res, err = quadgk(iq, -lim, lim)
    return res
end

function using_integrals(eta::Real, lim::Real)
    ip = IntegralProblem(integrand, -lim, lim, eta)
    return solve(ip, QuadGKJL())[1]
end

er = LinRange(-30, 2, 200)
lim = 100
plot(er, using_quadgk.(er, lim), yscale=:log10, label="QuadGK")
plot!(er, using_integrals.(er, lim), yscale=:log10, label="Integrals")

lim = 1000
plot(er, using_quadgk.(er, lim), yscale=:log10, label="QuadGK")
plot!(er, using_integrals.(er, lim), yscale=:log10, label="Integrals")

The function should be smooth, so the calculation from Integrals.jl doesn't make sense. I've also tested other solver algorithms, but they give me also wrong results in varying degree.

Julia 1.6.6 Integrals v3.0.2