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?
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.