Allow EnsembleProblem to re-use integrator cache

[PerezHz]

Is your feature request related to a problem? Please describe.

I'm doing many ODE integrations for different initial conditions in an EnsembleProblem with OrdinaryDiffEq. In this case the integrator cache is re-built every time for every iteration, causing allocations on every iteration, with more allocations for higher-order methods.

Describe the solution you’d like

It would be useful to be able to re-use integrator caches in the iteration of an EnsembleProblem instead of having to construct it every time.

Describe alternatives you’ve considered

As suggested by @ChrisRackauckas on Slack, a workaround is to use the integrator interface, i.e. doing __init then solve! and reinit!, instead of calling solve directly.

Additional context

A simplified example of the requested feature:

using OrdinaryDiffEq
f!(du, u, p, t) = (du .= u)
q0 = [0.2, 0.0, 0.0, 3.0]
prob = ODEProblem(f!, q0, (0.0, pi))
prob_func(prob, i, repeat) = remake(prob; u0 = q0 + 1e-8randn(4))
ensprob = EnsembleProblem(prob; prob_func)
sim = solve(ensprob, Vern9(); trajectories=10) # cache is re-built on every iteration