I need some help with this problem. I need to solve the differential equation $$xy''+3y'+x^3y=0$$ using power series. I used Frobenius method to expand about $x=0$ since it is an singular regular point. So I assumed a solution $y(x)=\sum_{j=0}^\infty a_jx^{s+j}$. After replacing it on the origina equation I ended up with: $$\sum_{j=0}^\infty a_j(s+j)(s+j-1)x^{s+j-1}+3\sum_{j=0}^\infty a_j(s+j)x^{s+j-1}+\sum_{j=0}^\infty a_jx^{s+j+3}=0$$ In order to equal the exponents of $x$, I expanded the first four terms in the first two series and ended up with the inditial equation: $$s^2+2s=0\Rightarrow s=0 \qquad s=-2$$ Now, because of the fact that the two roots differ by an integer, that means that the higher root will yield a solution while the smaller may or may not. Thus, for $s=0$ I found that the solution is $$y_1(x)=a_0\sum_{j=0}^\infty\frac{(-1)^j}{2^{2j}(2j+1)!}x^{4j}$$ I'm struggling to find the second solution. I tried to use two methods. The first one by using: $$y_2(x)=y_1(x)\int^x\frac{\exp\left[-\int^{x_2}P(x_1)dx_1\right]}{\left[y_1(x_2)\right]^2}dx_2$$ where $P(x)=\frac{3}{x}$, but I dont know what to do with the term $[y_1(x)]^2$.
The second method I tried was to use the series form of the second solution that my book (Mathematical Methods for physicists, Arfken) gives by writing $P(x)=\sum_{i=-1}^\infty p_ix^i$ and $Q(x)=\sum_{j=-2}^\infty q_jx^j$, replacing that in the integral form of the first method gives $$y_2(x)=y_1(x)\ln|x|+\sum_{j=-n}^\infty d_jx^{j+\alpha}$$ where $n$ is the difference between the two roots of the inditial equation and $\alpha$ is the higher root.
When I tried this method, I replaced $y_2(x)$ in the original ODE usiing $n=2$ and $\alpha=0$. After taking the derivatives and with the fact that $y_1(x)$ is a solution, I ended up with this: $$2\left(\frac{y_1(x)}{x}+y_1'(x)\right)+\sum_{j=-2}^\infty [j(j-1)+3j]d_j x^{j-1}+\sum_{j=-2}^\infty d_j x^{j+3}=0$$ The problem here is the term $2\left(\frac{y_1(x)}{x}+y_1'(x)\right)$, I don't know what to do with it. In my book, they write it as a new power series, but I don't know how to determine the coefficients of such power series.
I apologize for the long post, but I wanted to show a bit of the process I did, Hope you can help me. If it is necessary, I can upload an image of the whole process I did.
In your first approach, $y_2(x)=y_1(x)\int^x\frac1{s^3y_1(s)^2}ds$, you need to take into account that $y_1(x)$ is a series in $x^4$, so that the same is the case for $y_1(x)^{-2}=1+b_4x^4+b_8x^8+...$ Inserting that gives $$ \int^x\frac1{s^3y_1(s)^2}ds=\int^x(s^{-3}+b_4s+b_8s^5+...)ds=-\frac12s^{-2}+\frac12b_4s^2+\frac16b_8s^6+... $$ so that for this equation you do not get logarithmic terms, both basis solutions are Frobenius power series. This means that you can go back to the start and compute the second basis solution the same way as the first one via the coefficient recursion.