Solutions to an ODE using Frobenius when roots are repeated I have been looking at the following differential equation
$$xy^{\prime \prime} + y^{\prime} + xy=0$$
Observe that since $x=0$ is a regular singular point we use the Frobenius method and substitute $y=\sum_{n=0}^{\infty}a_n x^{n+r}$ to get
$$r^2a_0x^{r-1} + (r+1)^2a_1x^r + \sum_{n=2}^{\infty}[(n+r)^2a_n+a_{n-2}]x^{n+r-1}=0$$
Since $a_0 \neq 0$ we get $r=0,0.$ giving $y_1 = a_0(1-\dfrac{x^2}{2^2}+\dfrac {x^4}{4^2\cdot2^2}-\dfrac{x^6}{6^2\cdot4^2\cdot2^2}......)$.
Now apparently to obtain the second solution I have to differentiate this solution w.r.t. $r$ and then evaluate this at $r=0$, but I really don't see where this comes from. Could anyone show me explicitly why exactly we do this?  Or more generally -

If the roots of the indicial equation $r_1$ and $r_2$ are repeated,  why is the second solution of the form
$$y_2= c\left(\frac {\partial y_1}{\partial r}\right)_{\large {r=r_1}}?$$

All the sources I've looked at just state this fact without really showing where this came from, so either I'm missing something really obvious or there is some other explanation.
 A: If you modify the equation to
$$
xy''+(1-\varepsilon)y'+xy=0,~~~ε\approx 0
$$
then the coefficient equation changes to
$$
x^{n-1+r}:~~ (n+r)(n+r-ε)a_n+a_{n-2}=0
$$
The indicial equation at $n=0$ now has two different solutions $r=0$ and $r=ε$ that have no integer distance, so both give a Frobenius power series solution. $a_k=0$ for odd $k$ results in both cases. The coefficient recursions for the even-index coefficients, $n\ge 2$ even, are
$$
r=0:~~a_n=-\frac{1}{n(n-ε)}a_{n-2}
\\
r=ε:~~a_n=-\frac{1}{(n+ε)n}a_{n-2}
$$
Thus we get the first solution $y_1(x)=y(-ε,x)$ and the second $y_2(x)=x^εy(ε,x)$, where $y(ε,x)$ was set as the power series with the coefficients of the second solution with $a_0=1$. In the limit $ε\to 0$ both give the same solution.
Now as also linear combinations of both are solutions, we get
$$
y_3(x)=\frac{y_2(x)-y_1(x)}{ε}=\frac{x^ε-1}{ε}y(ε,x)+2\frac{y(ε,x)-y(-ε,x)}{2ε}
$$
as a third solution. In the limit the difference quotients become differential quotients, the limit is
$$
y_3(x)=\ln(x)y(0,x)+2\left.\frac{\partial y(ε,x)}{\partial ε}\right|_{ε=0}.
$$
This shows how a derivative can be introduced into the solution process. Now try to map this to the somehow parametrized coefficient sequence of the unmodified solution. I think it will work, but I see not how you could get a direct justification from that. In the direct interpretation of what you wrote, you would consider
$$
y(r,x)=x^r\sum_{n=0}^\infty a_n(r)x^k
$$
where
$$
a_n(r)=\frac1{(n+r)^2}a_{n-2}
$$
so that $y_1(x)=y(0,x)$ and the $r$-derivative of $y(r,x)$ at $r=0$ results in the same formula for $y_3$ at $r=ε=0$.
