I am trying to solve $$xy''+2xy'+6e^xy=0,\space x>0$$ about $x=0$ using Frobenius method. The question is from Boyce & DiPrima.
I have found the exponents at singularity to be $r_1=1$ and $r_2=0$. I have also found the series solution for $r=1$ by substituting $$y=x^r\sum_{n=0}^{\infty}a_nx^n$$ into the ODE and replacing $e^x$ with its Taylor series:
$$x^r\sum_{n=-1}^\infty a_{n+1}(n+r+1)(n+r)x^n+2x^r\sum_{n=0}^\infty a_{n}(n+r)x^n + 6x^r\sum_{n=0}^\infty (a_{n}x^n) \times \sum_{n=0}^\infty \frac{x^n}{n!} = 0$$ $$\dots$$ $$a_1 = -a_0\frac{2r+6}{r(r+1)}$$ $$a_2 = -a_0\frac{-\frac{(2r+8)(2r+6)}{r(r+1)}+6}{(r+2)(r+1)}$$ $$\dots$$
When I set $r=1, a_0=1$ and calculate these coefficients, I get $y_1=x\left(1-4x+\frac{17}{3}x^2-\frac{47}{12}x^3+\dots\right)$ which is correct (according to my book).
Now, to find the second solution (associated with $r_2=0$), I use this equation: $$y_2=c\ln(x)y_1(x) + x^{r_2}\sum{b_nx^n} \tag{1} \label{y2-form}$$
I find $$c = \lim_{r\rightarrow r_2}(r-r_2)a_1(r)=-\frac{2r_2+6}{r_2+1}=-6.$$
This is also correct. But my results for $b_n$ are wrong. To find them, I set $a_0=r-r_2=r$, find $a_n(a_0=r, r)$ (by canceling the factors $r$), and differentiate them with respect to $r$: $$b_n=\left.\frac{d a_n(a_0=r, r)}{d r}\right\vert_{r=0}$$
$$a_0 = r$$ $$b_0 = \left.\frac{d r}{d r}\right\vert_{r=0} = 1$$ $$b_1 = \left.\frac{d}{d r}\left(-\frac{2r+6}{r+1} \right)\right\vert_{r=0} = 4$$ $$b_2 = \dots = -49$$
These are completely wrong. The books says $b_0 = 1, b_1 = 0, b_2 = -33, b_3 = \frac{449}{6}x^3, \dots$
I have seen this method of finding $b_n$ in a proof of Frobenius method in the book An Introduction to Ordinary Differential Equations written by Coddington. Did I get all this "multiplying by $(r-r_2)$ and differentiating" thing wrong? Or maybe I am doing a silly mistake. Boyce & DiPrima suggests putting $\eqref{y2-form}$ in the ODE and calculating the coefficients that way. But I want to understand what is wrong with the other method.
Since $y_1$ and $y_2$ are linearly independent, we can add any multiple of $y_1$ to $y_2$, and get a solution that is still independent with $y_1$. Since $y_1$ is a power series starting with $x^1$, this means that the coefficient $b_1$ in the second solution is arbitrary.
Notice that the difference between the $b_n$ values found in the question and the $b_n$ values in the book are multiples of the coefficients in $y_1$. So we obtain the book's solution by subtracting $4y_1$ from $y_2$. This means that the $y_2$ found here is also correct.