On the wikipedia explanation of Frobenius method and the indical equation I cant manage to understand the "explanation" on wikipedia for the Frobenius method.(https://en.wikipedia.org/wiki/Frobenius_method)
Why does the solution one constructs satisfy
$z^2u_{r}''+tPu_{r}'+Qu_{r}=I(r)z^r$ ?
I can see why the first coefficient is zero when $I$ is but not the rest of them. Or maybe there is another explanation why the $I$ appears on the RHS.
 A: Well, it seems to me that the wikipedia article is raising much more questions than it is answering. It is quite sad that it does not reflect any of the more recent viewpoints on the Frobenius method instead of merely reproducing some old-fashioned ugly formulas. Nevertheless, if you want to gain a really detailed insight of what is happening here (although quite a bit messy still) in case of the second-order equations I strongly recommend the book by a certain Marvin Goldstein to be got for free from here https://eric.ed.gov/?id=ED089982 I shall attempt to answer your question from slightly different point of view, taking inspiration from the lecture notes http://people.math.umass.edu/~cattani/hypergeom_lectures.pdf (containing a lot of misprints, so that you should check everything for yourself before use) and from the book [IKSY]: https://books.google.cz/books?id=hCcGCAAAQBAJ&pg=PA328&lpg=PA328&dq=iwasaki+kimura&source=bl&ots=rghGvRV-HP&sig=UEyrYmBmvjjAnbG83Gl3ooNA-rs&hl=cs&sa=X&redir_esc=y#v=onepage&q=iwasaki%20kimura&f=false (the other well known standard references like coddington and levinson or ince should do as well for most of what we need)
So suppose we are given an $n-$th order ordinary differential equation
$$
Lw := (a_0(z)\partial^n + a_1(z)\partial^{n-1} + \ldots + a_{n-1} \partial + a_n(z))w = 0
$$
with the unknown function $w$ and $\partial$ to be the derivative with respect to the only independent variable, say $z.$ Here it is supposed that $a_j(z)$ are holomorphic in an open set $U \subset \mathbb{C}.$ Moreover, let $z_0 \in U$ be an isolated regular singular point of that equation meaning that although $a_0(z_0) = 0,$ the functions $a_j(z)/a_0(z)$ have at worst poles of order $j$ at $z_0$ meaning that after multiplying the equation
$$
\partial^n w + \frac{a_1(z)}{a_0(z)}\partial^{n-1} w + \ldots + \frac{a_{n-1}(z)}{a_0(z)} \partial w + \frac{a_n(z)}{a_0(z)} w = 0
$$
through by $z^n,$ you obtain the equation
$$
z^n\partial^n w + b_1(z)z^{n-1}\partial^{n-1} w + \ldots + b_{n-1}(z)z \partial w + b_n (z) w = 0,
$$
where $b_j$'s are holomorphic at $z_0.$ We can suppose without any loss of generality that $z_0 = 0.$ Now we introduce the operator $\theta = z \partial.$ This is a very handful shortcut. You can check for yourself that $\theta(z^{\alpha}) = \alpha z^{\alpha}$ for any complex $\alpha$ (with a suitable choice of the branch of the argument) and that $z^m \partial^m =\theta (\theta - 1) \ldots (\theta - m +1),$ so that in terms of this $\theta$ the last ODE reads
\begin{eqnarray}
&&\theta (\theta - 1) \ldots (\theta - n +1)w + b_1 (z) \theta (\theta - 1) \ldots (\theta - n +2)w + \ldots + b_{n-1}(z) \theta + b_n(z) \\
&=& \sum_{i=0}^{n} b_{n-i}(z) \theta (\theta - 1) \ldots (\theta - i +1) = 0,
\end{eqnarray}
where of course $b_0(z) = 1.$ Now, since the $b_{n-i}$'s are still holomorphic at zero, we have the Taylor expansions
$$
b_{n-i}(z) = \sum_{j=0}^{\infty}\frac{b_{n-i}^{(j)}(0)}{j!} z^j
$$
for every $i=0,\ldots, n.$ Seek a solution in the form $w(z) = z^{\lambda} \sum_{k=0}^{\infty} u_k z^k.$ Plugging everything into our $Lw = 0$ (more precisely into its most recent form with $\theta$s) and manipulating series a bit (and using $\theta z^j = jz^j$ many many times) we obtain:
\begin{eqnarray}
&&\sum_{i=0}^{n}\sum_{j=0}^{\infty} \frac{b_{n-i}^{(j)}(0)}{j!} z^j \theta (\theta - 1) \ldots (\theta - i +1) \sum_{k=0}^{\infty} u_k z^k \\
& = & \sum_{k=0}^{\infty} \sum_{j=0}^{\infty} \sum_{i=0}^{n} u_k \frac{b_{n-i}^{(j)}(0)}{j!} z^j \theta (\theta - 1) \ldots (\theta - i +1) z^{\lambda + k} \\
& = & \sum_{k=0}^{\infty} \sum_{j=0}^{\infty} u_k z^j \left( \sum_{i=0}^{n} \frac{b_{n-i}^{(j)}(0)}{j!} \theta (\theta - 1) \ldots (\theta - i +1) \right) z^{\lambda + k}\\
& = & \sum_{k=0}^{\infty} \sum_{j=0}^{\infty} u_k z^j \left( \sum_{i=0}^{n} \frac{b_{n-i}^{(j)}(0)}{j!} (\lambda+k) (\lambda + k - 1) \ldots (\lambda + k - i +1) \right) z^{\lambda + k}.
\end{eqnarray}
Now we define
$$
p^{(j)}(\lambda):= \sum_{i=0}^{n} \frac{b_{n-i}^{(j)}(0)}{j!} \lambda(\lambda - 1) \ldots (\lambda - i +1).
$$
Note that for $j=0$ we have the usual indicial polynomial $p(\lambda).$ This means that the penultimate equation above can be written as
\begin{eqnarray}
&& \sum_{k=0}^{\infty} \sum_{j=0}^{\infty} u_k z^j \left( \sum_{i=0}^{n} \frac{b_{n-i}^{(j)}(0)}{j!} (\lambda+k) (\lambda + k - 1) \ldots (\lambda + k - i +1) \right) z^{\lambda + k} \\
& = & \sum_{k=0}^{\infty} \sum_{j=0}^{\infty} u_k z^j p^{(j)}(\lambda + k) z^{\lambda + k} \\
& = & \sum_{k=0}^{\infty} \sum_{j=0}^{\infty} u_k z^{k+j} p^{(j)}(\lambda + k) z^{\lambda}.
\end{eqnarray}
Making the change of variables (for this technique I refer to the splendid https://books.google.cz/books/about/Special_Functions.html?id=z-ApAQAAMAAJ&redir_esc=y and its Chapter 4, section 37) $k+j \rightarrow k,$ meaning that you define new summing variables, say $l$ and $m,$ connected with the old ones by the formula $j = l$ and $k = m-l,$ which has the effect of shifting the domain of summation to $m\geq 0$ and $0 \leq l \leq m$ and then you write back $k$ instead of $m$ and $j$ instead of $l,$ the variables being dummy, you obtain
\begin{equation}
\sum_{k=0}^{\infty} \sum_{j=0}^{\infty} u_k z^{k+j} p^{(j)}(\lambda + k) z^{\lambda} = \sum_{k=0}^{\infty} \sum_{j=0}^{k} u_{k-j}p^{(j)}(\lambda + k-j) z^{\lambda + k}.
\end{equation}
This means that the coefficients standing at $z^{\lambda + k}$ are given by
\begin{eqnarray}
&&\sum_{j=0}^{k} u_{k-j}p^{(j)}(\lambda + k-j) = u_k p(\lambda + k) + \sum_{j=1}^{k} u_{k-j}p^{(j)}(\lambda + k - j) \\
& = & u_k(\lambda + k) + \sum_{j=0}^{k-1} u_{k-j+1}p^{(j+1)}(\lambda + k - j+1)
\end{eqnarray}
and the coefficient $u_k$ with $k = 0$ has to satisfy
$$
u_0p(\lambda) = 0
$$
if the sought function $w$ is to be at least a formal solution to the ODE above (note that we are not talking about convergence yet!). The point now is that if $p(\lambda + k) \neq 0$ for every $k \in \mathbb{N}$ and $u_0 \neq 0$ (say $u_0=1$) then we can find the $u_k$'s recursively starting with $u_0$ in such a way that the expression
$$
\sum_{j=0}^{k} u_{k-j}p^{(j)}(\lambda + k - j)
$$
vanishes for every $k \geq 1$ and this is clearly just for this choice of the coefficients $u_k$ that we have
$$
L(w) = p(\lambda) z^{\lambda}
$$
(this should be pretty obvious now). Taking the special case $n=2$ back again could perhaps make things clearer.
