Incorrect results with TMJets for initial condition containing an equilibrium point

[schillic]

(copied here from another discussion)

using ReachabilityAnalysis, DifferentialEquations, Plots

@taylorize function f!(dy, y, p, t)
    dy[1] = y[1] * cos(y[1])
end

prob = @ivp(y' = f!(y), y(0) ∈ 0 .. 2, dim=1)

alg = TMJets21b(orderT=10, orderQ=2, abstol=1e-15)

sol = solve(prob, tspan=(0.0, 5.0), alg=alg,
            ensemble=true, trajectories=200, adaptive=true, dt=0.001,
            trajectories_alg=Rodas5())

# plot flowpipe and the ensemble solution
plot(sol, vars=(0, 1), linewidth=0.2, xlab="t", ylab="y(t)")
plot!(ensemble(sol), vars=(0, 1), linealpha=1.0)

Something odd happens near 0. Mincing the initial interval gives answers that look better

prob = @ivp(y' = f!(y), y(0) ∈ split(Interval(0, 2), 50), dim=1)
alg = TMJets21b(orderT=10, orderQ=2, abstol=1e-15)
solm = solve(prob, tspan=(0.0, 5.0), alg=alg)

# plot flowpipe and the ensemble solution
plot(solm, vars=(0, 1), linewidth=0.2, xlab="t", ylab="y(t)")
plot!(ensemble(sol), vars=(0, 1), linealpha=1.0)

If the interval of interest is 1 .. 2, or any other not containing 0, the integrator behaves properly.

prob = @ivp(y' = f!(y), y(0) ∈ Interval(1, 2), dim=1)
alg = TMJets21b(orderT=10, orderQ=2, abstol=1e-15)
solm = solve(prob, tspan=(0.0, 5.0), alg=alg)

# plot flowpipe and the ensemble solution
plot(solm, vars=(0, 1), linewidth=0.2, xlab="t", ylab="y(t)")
plot!(ensemble(sol), vars=(0, 1), linealpha=1.0)

There is something about the initial interval having an (unstable) equilibrium point inside. If we include 0 in the initial interval, at some point things go wrong. If 0..0 is the initial point for the expansion, so the initial interval is symmetric around zero, then we get an error: "Minimum tolerance reached" before finishing the integration. If the interval is asymmetric, the integration may finish but the flowpipe obtained is wrong, or simply the error bounds of the Taylor expansion become useless (-Inf .. Inf).