Large error when integrating around point of discontinuity

Describe the bug 🐞

The numerical solution has a large discrepancy from the analytical solution when the bound is close to the point of discontinuity.

Expected behavior

I expect the difference should be really small (much smaller than 0.00059).

Minimal Reproducible Example 👇

using Integrals

function f(x)
    if x < 1
        return exp(x)
    elseif x >=1
        return exp(2x)
    end
end

function F(x0, x1)
    if x0 < 1 && x1 < 1
        return exp(x1) - exp(x0)
    elseif x0 < 1 && x1 >= 1
        return exp(2x1) / 2 - exp(2) / 2 + exp(1) - exp(x0)
    elseif x0 >= 1 && x1 >= 1
        return exp(2x1) / 2 - exp(2x0) / 2
    end
end


x0 = 0.997
x1 = 2.0

prob = IntegralProblem((x,p) -> f(x), (x0,x1))
sol1 = solve(prob, QuadGKJL(), reltol=1e-12, abstol=1e-12).u
sol0 = F(x0, x1)

sol1 - sol0
(sol1 - sol0) ./ sol0

Environment (please complete the following information):

• Output of using Pkg; Pkg.status()

[de52edbc] Integrals v4.6.0

• Output of versioninfo()

Julia Version 1.11.5
Commit 760b2e5b739 (2025-04-14 06:53 UTC)
Build Info:
  Official https://julialang.org/ release
Platform Info:
  OS: macOS (arm64-apple-darwin24.0.0)
  CPU: 14 × Apple M4 Pro
  WORD_SIZE: 64
  LLVM: libLLVM-16.0.6 (ORCJIT, apple-m1)
Threads: 1 default, 0 interactive, 1 GC (on 10 virtual cores)

Reproduce the minimal example using IntegralProblem and QuadGKJL() with the stated tolerances, then inspect the integration path for bounds near the discontinuity at x=1. Compare the numerical result with F(x0, x1); done means the discrepancy is substantially smaller than 0.00059 or the issue is clearly explained and covered by a regression test.