Large interpolation error for exponential integrators in DifferentialEquations

I was recently working with `SplitODEProblem` using exponential integrators and I noticed a huge error when doing the interpolation of the solution, which I think is wrong.

I began to notice this since, by default, the solution gives the values at each time step and always ends at `tspan[2]`, regardless if this value is smaller than the last step for a given `dt`. This last interpolated step seemed to be way off from the solution, so I started experimenting with it and I have the following MWE.

``````using DifferentialEquations, GLMakie

K = [1.0 -1.0 0.0;
-1.0 2.0 -1.0;
0.0 -1.0 2.0]
M = [1.0 0.0 0.0;
0.0 1.0 0.0;
0.0 0.0 1.0]

F = [0.1, 0.0, 0.0]
F₀ = [20.0, 20.0, 0.0]

M⁻¹F, M⁻¹F₀, M⁻¹K = M \ F, M \ F₀, M \ K

ξ = 0.01

function f2!(du, u, p, t)
N = length(u) ÷ 2
M⁻¹F, M⁻¹F₀, Ω = p

du[(N+1):end] .= M⁻¹F * sin(Ω * t) .+ M⁻¹F₀
end

A = MatrixOperator([zeros(size(M)) I; -M⁻¹K -M⁻¹K*ξ])

Ω = 0.445
p = (M⁻¹F=M⁻¹F, M⁻¹F₀=M⁻¹F₀, Ω=Ω)

tspan = (0.0, 2π / Ω * 10)
u0 = zeros(6)
prob = SplitODEProblem(A, f2!, u0, tspan, p)
sol = solve(prob, LawsonEuler(), dt=1e-1, progress=true)

t = (2π/Ω*9):(2π/Ω/100):(2π/Ω*10)
lines(t, sol(t)[1, :])
``````

This example creates an oscillating dynamical system and I solve it with an exponential integrator. However, if I do the interpolation as in the last lines, I get the following plot

whereas examining the complete solution from the integrator gives

which seems to be smooth and well resolved. Therefore, is this a bug of the interpolation functions for these kinds of scheme? Also, is it intended that the last time step from the time vector is always `tspan[2]` even though there could be a time instant that goes after this point?

Can you open an issue?

Done, you can find it here.