Does the Frobenius method work for all second order linear differential equations with only regular singular points? For $$x^2y''+x(x^2+1)y'+(x-4)y=0,\tag1$$ there is a regular singular point at $0$, but when I tried to use the Frobenius method and substituted $$y=\sum_{n=0}^\infty a_n x^\left(n+r\right)\tag2$$ into the differential equation, I obtained 
\begin{multline}
(r^2-4)a_0x^r+[(r^2+2r+5)a_1+a_0]x^\left(r+1\right)
\\+\sum_{n=0}^\infty\left\{[(n+r+2)^2+2]a_\left(n+2\right)+a_\left(n+1\right)+(n+2)a_n\right\}=0\tag3
\end{multline} 
I understand that if all $x$ in an interval satisfy the differential equation, then all the coefficients must be zero. But then this means $r=\pm2$ (from the first term), and certainly the second term will not be satisfied by arbitrary $a_0$ and $a_1$ (the initial conditions). 
Therefore, my question is does the Frobenius method only work for certain second order linear differential equations with only regular singular points, like where $p(x)$ and $q(x)$ in $~x^2y''+p(x)y'+q(x)y=0~$ are first or second degree polynomials?
 A: You are selecting the values of $r$ specifically so that the first term is zero for any arbitrary value of $a_0$.
You made some sign errors in the collection of coefficients. Using $a_{-2}=a_{-1}=0$, the insertion of the series expansion can be written shortly as
$$
\sum_{n=0}^\infty\left\{[(n+r)^2-4]a_n+(n+r-2)a_{n-2}+a_{n-1}\right\}x^{n+r}=0\tag1
$$
Then the first two coefficient identities are
$$
[r^2-4]a_0=0\\
[(r+1)^2-4]a_1+a_0=0\tag2
$$
so that still $r=\pm2$, but now  $a_1=-\frac1{2r+1}a_0$. Now the question remains if all other coefficients are determined by $a_0$. For $r=2$ this is certainly the case, but for $r=-2$ one gets for $n=1,..,4$ the equation
$$
-3a_1+a_0=0\\
-4a_2+a_1-2a_0=0\\
-3a_3+a_2-a_1=0\\
0a_4+a_3+0a_{2}=0\tag3
$$
which enforces $a_0=0$ so that this branch gives the same solution set as the case $r=2$.
You obtain the second solution by order reduction, as the first solution ($r=2$) has the form $y(x)=x^2a(x)$, $a(x)=\sum_{n=0}^\infty a_nx^n$, $a_0=1$, setting per the reduction-of-order method $y_2=u(x)y(x)=x^2u(x)a(x)$ gives for $u$ the reduced equation
\begin{align}
y_2'(x)&=u'(x)y(x)+u(x)y'(x)\\
y_2''(x)&=u''(x)y(x)+2u'(x)y'(x)+u(x)y''(x)\\[1em]
\hline
0&=x^2[u''(x)y(x)+2u'(x)y'(x)]+x(x^2+1)u'(x)y(x)\tag4\\[1em]
\frac{u''(x)}{u'(x)}&=-\frac{2x^2y'(x)+x(x^2+1)y(x)}{x^2y(x)}=-\frac4x-\frac{2a'(x)}{a(x)}-x-\frac1x\tag5
\end{align}
so that with the most simple integration constant one finds
$$
u'(x)=\frac{e^{-x^2/2}}{x^5a(x)^2}\tag6
$$
While one could directly compute this per power series division and term-wise integration, the structure information alone tells us that
$$
u(x)=b_{-4}x^{-4}+...+b^{-1}x^{-1}+\ln(x)+b_1x+...\tag7
$$
so that
$$
y_2(x)=x^2\ln(x)a(x)+\sum_{n=0}^\infty c_nx^{n-2}.\tag8
$$

As to the general question, yes, at a regular singularity you get at least one solution in the form of a generalized power series, provided the roots of the indicial equation are real. If they are complex, then with $x^{a+ib}=x^a(\cos(b\ln x)+i\sin(b\ln x))$ one gets some more complicated terms in the real forms of the solution.
