Let $x\in R^n$ and consider the system
$$ \dot{x}=f(t,x) \;\;\mbox{with}\;\; x(0)=x_0 $$ and suppose that we know it's exact or very accurate solution $x(t)$ for the time interval $[0,T]$.
I'm asking myself how to approximate the value of a second solution at time $T$, say $y(T)$, of the same system for small variations of the initial conditions, i.e.
$$ \dot{y}=f(t,y) \;\;\mbox{with}\;\; y(0)=x_0 + \epsilon v $$
for do that I think that it is necessary to compute the derivative of the state w.r.t. the initial conditions, so, lets write each component of $x$ as $x^i$, the same for $f$.
$$ \frac{\partial }{\partial x_0^j} \dot{x}^i = \frac{\partial f^i}{\partial x^k} \frac{\partial x^k}{\partial x_0^j} $$
for shorthand, I call $\psi^i_j = \frac{\partial x^i}{\partial x_0^j}$. So the solution of the system $$ \dot{\psi}_j^i = \left[\frac{\partial f^i}{\partial x^k}\right]_{x=x(t)} \psi^k_j $$ and the equation $$ y^i(t) = x^i(t) + \psi^i_j v^j $$ should give the answer to my question for all $t$. Is that true?, second, what should be the initial conditions for $\psi_j^i$?, and third, it is that procedure correct for the variation of the solution w.r.t. other parameters if they depends on time?