DifferentialEquations.jl solution object: save_idxs and interpolation

using DifferentialEquations
N = 10

H = randn(N,N)
H -= H'

ψ0 = zeros(N); ψ0[1] = 1

f(ψ :: Array{Float64,1},p,t) = H*ψ

prob = ODEProblem(f, ψ0, (0.0,1.0), save_idxs=[1])
sol = solve(prob, Tsit5(), reltol=1e-8, abstol=1e-8)
sol(0.5)

is fine. Or using saveat is fine. Vern7(lazy=false) would probably be a good choice here. It’s an obscure issue that comes from the exact combination of features you chose and the algorithm that got automatically picked. It’s kind of amazing you found it :wink:. We’ll try to make this into account in the default algorithm chooser.