The Effect of Small Variations in the Coefficients of a Linear Differential System In E.L. Ince's Ordinary Differential Equations, page $219$, the effect of small variations in the coefficients of a linear differential system is examined. In particular, the main result of the section is that the "index" of the system is not raised by any variation of the coefficients which is "uniformly small."
The result above concerns the "index" of the system, but in general, the question of the effect of small variations in the coefficients of a differential equation sounds interesting and, despite going over a couple of texts, modern and old, I don't see anyone discussing about it. I have one specific question in mind: Suppose we focus on just ordinary linear differential equations and one (or maybe all) of the terms is differentiably parametrized (i.e. the partial derivative of the term w.r.t. the parameter exists - or any other way you want to interpret it). Is the solution (now a two-variable function) differentiable w.r.t. the parameter as well? Maybe the answer to this question can be found in PDE's, but I hardly know any PDEs sensibly.
In general, does anyone know of a good text (preferably a modern one) that addresses such a question?
Thanks!
Note: By the way, I am not sure if this kind of question is allowed here, but judging from how this topic evades internet searches, if anyone is interested to look into the problem deeper and from other angles (e.g. computational, applied), I would like to do the work together - preferably under a formal position such as a postdoc.
 A: I don't know anything about relevant literature but here are some things you can say with your bare hands. Consider the second-order homogeneous linear ODE
$$\frac{d^2 y}{dt^2} + b \frac{dy}{dt} + c y = 0$$
with fixed initial conditions $y(0) = 0, y'(0) = 1$. We know, very classically, that the solutions to this ODE behave drastically differently depending on the behavior of the roots of the characteristic polynomial $\lambda^2 + b \lambda + c$; let's first compute the solution with fixed initial conditions and then see what happens if we vary $b$ and $c$. Write
$$r_1 = \frac{-b + \sqrt{b^2 - 4c}}{2}, r_2 = \frac{-b - \sqrt{b^2 - 4c}}{2}$$
for the two roots of the characteristic polynomial. I will only write down the solution in the case $r_1 \neq r_2$ because it turns out a limit will get us the repeated root case. The solution is a linear combination
$$y(t) = A_1 \exp(r_1 t) + A_2 \exp(r_2 t)$$
and checking compatibility with initial conditions gives $A_1 + A_2 = 0, A_1 r_1 + A_2 r_2 = 1$, which gives
$$\boxed{ y(t) = \frac{\exp (r_1 t) - \exp (r_2 t)}{r_1 - r_2} }.$$
This actually gets us the repeated root $r_1 = r_2$ case: we can take the limit as $r_2 \to r_1$ and get
$$y(t) = t \exp (rt)$$
where $r = r_1 = r_2$ is the common repeated root. This is a cute way to understand where the extra factor of $t$ comes from.
In this example $y(t)$ depends smoothly as a function on the parameters $r_1, r_2$ but does not depend smoothly as a function on the original coefficients $b, c$, owing to the square-root terms in the quadratic formula; to see this more clearly we can set $b = 0$ so that $r_1, r_2 = \pm \sqrt{-c}$, giving
$$y(t) = \begin{cases} \frac{\exp(\sqrt{-c} t) - \exp(-\sqrt{-c} t)}{2 \sqrt{-c}} & \text{ if } & c \neq 0 \\ t & \text{ if } & c = 0 \end{cases}$$
We can make the answer look nicer by writing it in terms of $d = |\sqrt{c}|$: then we just have $y(t) = \frac{\sin dt}{d}$ if $c > 0$, $y(t) = t$ if $c = 0$, and $y(t) = \frac{\sinh dt}{d}$ if $c < 0$. Here are some plots of these functions for $d = 0, \pm \frac{1}{2}, \pm 1$ so you can get a sense of how they deform into each other, from WolframAlpha:

