Possibly generalize the indices?

[ChrisRackauckas]

I see lots of codes like:

for k in 1:sol.nd

Note that this will only loop though Array in Julia. If someone uses an OffsetArray or many other things in AbstractArray, this will not get the indices right, which is why it's normally nice to use eachindex(u) where u is the array if you're updating each element.

Do you plan on supporting non-traditional arrays? Note that this is something we'll be compiling a compatibility list of, so it's definitely okay to not go out of your way to do this. I think only OrdinaryDiffEq.jl is compatible with most of the odd types like this, while things like Sundials.jl, LSODA.jl, and ODEInterface.jl are on the other end and need Vector{Float64}.

[michakraus]

I have been thinking about this and it would be nice to generalise. However, I do not yet see an application for AbstractArrays other than Vector and Matrix, so the activation energy for me to look into this has not yet been high enough to cross the potential well.

[ChrisRackauckas]

I think some of the ODE setups from this thread are good motivators:

https://discourse.julialang.org/t/get-values-of-symbols-inside-a-function

or those which are multi-scale models

https://github.com/JuliaDiffEq/MultiScaleModels.jl

Full support of large swath of AbstractArrays is very time consuming though. I must admit I spent quite a bit of time to make it work without major performance losses, since there are many places you have to substitute out copy for deepcopy, save full type information, perform type checks, etc. While I think it ends up being too much to really commit to all at once, switching to using eachindex to do general indexing is a good first step and will already make OffsetArrays supported (provided you build other arrays correctly using zeros or similar and passing the indices, so I guess there is more too it...)