Optimizing the computation of jacobians for multiple shooting implementation

[gdalle]

@jackyarndley

I am primarily wondering if there is a way to avoid the repeated cost of __init calls when solve is called in the context of using AD on the solver, as using the reinit! function tended to be slower.

That might be more of a question for BoundaryValueDiffEq.jl developers then