Using `saveat` with IDA results in derivatives with sharp discontinuities

[topolarity]

Describe the bug 🐞

If saveat is provided to solve() when using IDA, the resulting derivative values appear to have strong discontinuities.

Minimal Reproducible Example 👇

julia> using OrdinaryDiffEq, Sundials, Plots
julia> dae_f(du, u, _, t) = [du[1] - u[2], u[2] - sin(t)]
julia> prob = DAEProblem(DAEFunction(dae_f), [0.0, 0.0], [0.0, 0.0], (0., 2π), nothing,
                         differential_vars=[true,false], abstol=1e-6, reltol=1e-1)
julia> sol = solve(prob, IDA(), saveat=.01)
julia> plot(sol(sol.t, Val{1}))

Environment (please complete the following information):

• Output of using Pkg; Pkg.status()

  [459566f4] DiffEqCallbacks v2.36.1
⌃ [1dea7af3] OrdinaryDiffEq v6.70.0
  [c3572dad] Sundials v4.23.2

Additional context

Kudos to @oscardssmith for finding this one.

[topolarity]

Without saveat the interpolation is also quite cursed:

julia> sol = solve(prob, IDA())
julia> t′ = 0.0:1e-3:2π
julia> plot(sol(t′, Val{1}))

[ChrisRackauckas]

This seems pretty inherent to Sundials since we don't (can't?) recreate its internal polynomial.