Suppose $\mathbf{u}^1\in\mathbb{R}^n$ is an unknown and $\mu^1\in\mathbb{R}$ is a known parameter. Suppose we have solved the non-linear system of equation for $\mathbf{u}^1$ using Newton's iterations.
$$ \mathcal{N}(\mathbf{u}^1;\mu^1)=0 $$
Using $\mathbf{u}^1$ as an initial guess for the Newton iterations, it is possible to accelerate the solution of the system
$$ \mathcal{N}(\mathbf{u}^2;\mu^2) = 0 $$
where $\mu^1$ and $\mu^2$ are sufficiently close. Now, is it possible to extend this to ODEs? For instance, if I have solved
$$ \frac{d\mathbf{u}^1}{dt} = \mathcal{M}(\mathbf{u}^1;\mu^1) $$
Is it possible to accelerate the solution of
$$ \frac{d\mathbf{u}^2}{dt} = \mathcal{M}(\mathbf{u}^2;\mu^2) $$
where $\mu^1$ and $\mu^2$ again are sufficiently close.
The only thing I can think of is to integrate both ODEs implicitly with the same time-step. Solve the first ODE, and use the solution at each time-step to accelerate Newton's iterations for the second. Can anything else be done?