Solution of 2nd order linear ODE with regular singular points, and complex exponents at singularity The steady state temperature distribution of a rod given by:
    \begin{equation}
        \frac{\textrm{d}p(x)y'}{\textrm{d}x} - y = 0,\; 0 \leq x \leq 1,\; \text{and} \;y(0) = 0,
    \end{equation}
where $y(x)$ is the steady state temperature distribution, and $p(x) = x^s$ is the spatially dependent conductivity of the rod, for some $s \in [0, 1)$. For $0 \leq s < 1$ what is the form of the solution close to $x = 0$ (i.e. the first the term of the series solution)? Does a solution exist for $s = 1$? $s > 1$?

First, simplify the equation into a recognizable form:
    \begin{align*}
        &\quad \frac{\textrm{d}p(x)y'}{\textrm{d}x} - y = 0 \\
        &\equiv p(x)y'' + p'(x)y' - y = 0 \\
        &\equiv y'' + \frac{p'(x)}{p(x)}y' - \frac{1}{p(x)} y = 0 
    \end{align*}
Since $p(x) = x^s$, thus $p'(x) = sx^{s-1}$ and $p''(x) = s(s-1)x^{s-2}$. Substituting:
    \begin{align*}
        &\quad y'' + \frac{p'(x)}{p(x)}y' - \frac{1}{p(x)} y = 0 \\
        &\equiv y'' + \frac{s}{x}y - \frac{s^2 - s}{x^2}y = 0
    \end{align*}
We can immediately recognize that $x = 0$ is a regular singular point since:
    \begin{align*}
        &\lim_{x\rightarrow 0}  \frac{sx}{x} = s
        &\lim_{x\rightarrow 0} -\frac{(s^2 - s)x^2}{x^2} = s - s^2
    \end{align*}
It is worth noting that these limits exist for all $s$. 
So, we look for solutions of the form:
    \begin{align*}
        y = \sum_{n=0}^{\infty} a_nx^{n+r}
    \end{align*}
where $r(r-1) + sr + (s - s^2) = 0 = r^2 + (s-1)r + (s - s^2)$ is the corresponding indical equation. Solving for $r$, we obtain:
    \begin{align}
        \nonumber &\quad r = \frac{-(s-1) \pm \sqrt{(s-1)^2 - 4(s - s^2)}}{2} \\
        \nonumber &\equiv r = \frac{-(s-1) \pm \sqrt{s^2 - 2s + 1 - 4s + 4s^2}}{2} \\
        &\equiv r = \frac{-(s-1) \pm \sqrt{5s^2 - 6s + 1}}{2}
    \end{align}
For what values of $s$ do we get imaginary $r$? That would be when the discriminant of the quadratic formula is less than zero: 
    \begin{align*}
        5s^2 - 6s + 1 < 0
    \end{align*}
Consider when $s = 0$:
    \begin{align*}
        &\quad 5s^2 - 6s + 1 = 0 \\
        &\equiv s = \frac{6 \pm \sqrt{36 - 20}}{10} \\
        &\equiv s = \frac{6 \pm \sqrt{16}}{10} \\
        &\equiv s = \frac{6 \pm 4}{10} \\
        &\equiv s = 1 \vee s = \frac{1}{5}\\
    \end{align*}
Since for $(s=0 \wedge d = 5s^2 - 6s + 1) \implies d = 1$, then we know that the discriminant is positive for all values of $s \leq \frac{1}{5}$ and $s \geq 1$, and we have:
    \begin{align*}
        r_1 &= \frac{-(s-1) + \sqrt{(s - 1)(s - \frac{1}{5})}}{2} \\
        r_2 &= \frac{-(s-1) - \sqrt{(s + 1)(s - \frac{1}{5})}}{2}
    \end{align*}
The general solution is of the form:
    \begin{equation*}
        y = c_1\sum_{n = 0}^{\infty} a_nx^{n + r_1} + c_2\sum_{n = 0}^{\infty} a_nx^{n + r_2}
    \end{equation*}
where $a_n$ might be complex if $\frac{1}{5} s < 1$
We know that $y(0) = 0$, therefore:
    \begin{align*}
        &\quad y(0) = 0 = c_1\sum_{n = 0}^{\infty} a_n0^{n + r_1} + c_2\sum_{n = 0}^{\infty} a_n0^{n + r_2} \\
        &\equiv 0 = c_1\sum_{n = 0}^{\infty} a_n + c_2\sum_{n = 0}^{\infty} a_n \\
        &\equiv 0 = (c_1 + c_2) \\
        &\equiv c1 = -c_2
    \end{align*}
The first term of this series is:
    \begin{align*}
        \nonumber &\quad c_1\left(a_0x^{r_1} - a_0x^{r_2}\right) \\
        \nonumber &\equiv c_1a_0\left(x^{\frac{-(s-1) + \sqrt{(s - 1)(s - \frac{1}{5})}}{2}} - x^{\frac{-(s-1) - \sqrt{(s - 1)(s - \frac{1}{5})}}{2}}\right) \\
    \end{align*}
Let $\frac{1 - s}{2} = \alpha$ and $\frac{\sqrt{(s - 1)(s - \frac{1}{5})}}{2} = \beta$
If $0 \leq s \leq \frac{1}{5}$:
    \begin{equation*}
        y \approx c_1a_0\left(x^{\alpha + \beta} - x^{\alpha - \beta}\right)
    \end{equation*}
If $\frac{1}{5} < s < 1$:
    \begin{equation*}
        y \approx c_1a_0\left(x^{\alpha + i\beta} - x^{\alpha - i\beta}\right)
    \end{equation*}
Recall that:
    \begin{align*}
        &\quad x^r = e^{rln(x)} \\
        &\implies x^{\lambda + i\mu} = e^{\lambda\ln(x)}e^{i\mu\ln(x)} \\
        &\equiv x^{\lambda + i\mu} = e^{\lambda\ln(x)}(\cos(\mu\ln(x)) + i\sin(\mu\ln(x))) 
    \end{align*}
So we have oscillatory solutions.
Thus: 
    \begin{equation*}
        y \approx c_1a_0x^{\alpha}\left(2i\sin(\beta\ln(x))\right)
    \end{equation*}
Since the limits determining if $x = 0$ is a regular point exist for all $s$, we can write a (non-trivial) Frobenius series solution for all values of $s$.

Questions:
1) I have not used information about $0 \leq s < 1$, which is bothersome. Did I miss an opportunity to use it to simplify the form of the first term?
2) Is my reasoning for why solutions exist for $s = 1$ and $s > 1$ sound?
 A: $$y''+\frac{s}{x}y'-x^{-s}y=0$$
This is an ODE of the Bessel kind, but not on standard form. In order to make it standard, let $y(x)=x^aF(b\:x^c)=x^aF(X)$. The goal is to transform it to :
$$F''(X)+\frac{1}{X}F'(X)-\left(1+\frac{\nu^2}{X^2}\right)F(X)=0$$
The calculus is arduous. The result is :
$$a=\frac{1-s}{2} \: ; \: b=\frac{2}{2-s} \: ; \: c=\frac{2-s}{2} \: ; \: \nu=\frac{s-1}{s-2}$$
The general solution of the Bessel function is :
$$F(X)=c_1 I_\nu(X)+c_2 K_\nu(X)$$
$I_\nu(X)$ and $K_\nu(X)$ are the mofified Bessel functions of first and second kind. So, the general solution is :
$$y(x)=C_1 x^{\frac{1-s}{2}}I_{\frac{s-1}{s-2}}\left(\frac{2}{2-s}x^{\frac{2-s}{2}}\right) + C_2 x^{\frac{1-s}{2}}K_{\frac{s-1}{s-2}}\left(\frac{2}{2-s}x^{\frac{2-s}{2}}\right)$$
The formula is valid any $s$, but some particular cases arise sometimes with simplification of the Bessel functions to functions of lower level :
Case $s=1$ then $y=C_1 I_0\left(2 x^{1/2}\right)+C_2 K_0\left(2 x^{1/2}\right)$
Case $s=0$ then $y=C_1 \sinh(x) +C_2 \cosh(x)$
Case $0<s<1$ then :
In case of $0\leq x\leq 1$ , the condition $y(0)=0$ implies $C_2=0$
$$y(x)=C x^{\frac{1-s}{2}}I_{\frac{s-1}{s-2}}\left(\frac{2}{2-s}x^{\frac{2-s}{2}}\right)$$
The first term of the series expansion around $x=0$ of the Bessel function is :
$$I_\nu(X)\simeq\frac{1}{2^\nu \Gamma(\nu+1)}X^\nu$$
$$y(x)\simeq C \frac{1}{(2-s)^{\frac{s-1}{s-2}} \Gamma(1+\frac{1-s}{2-s})}x^{1-s} $$
This is consistent with a direct approach : 
Suppose that the first term be on the form $y\simeq cx^p$ then $y'\simeq cpx^{p-1}$ and $y''\simeq cp(p-1)x^{p-2}$
Puting them into the ODE : $cp(p-1)x^{p-2}+scpx^{p-2}-cx^{-s+p}\simeq 0$
$p(p-1+s)-x^{2-s}\simeq 0$
This implies $p-1+s=0$ hense $p=1-s$ and $y\simeq cx^{1-s}$ already obtained above.
