# Finding a solution to this equation using Frobenius method

by climbon
Tags: equation, frobenius, method, solution
 P: 1,666 Using the standard convention of letting $a_0=1$, I get: $$a_1=0$$ $$a_2=-\frac{h a_0}{(2+c)(1+c)+1+m^2}$$ with the remaining odd a_n=0 and even a_n equal to: $$a_n=-\frac{a_{n-2}+a_{n-4}}{(n+c)(n+c+1)+1+m^2}$$ with a similar relation for the other root expressed in b_n and so I can write the solution as: $$y(x)=K_1 \sum_{n=0}^{\infty} a_n x^{n+mi}+K_2 \sum_{n=0}^{\infty} b_n x^{n-mi},\quad x>0$$ and I believe because of the conjugates in the solution, for real initial conditions $y(x_0)=y_0$, and $y'(x_0)=y_1$, the imaginary component is annihilated leaving the desired real solution. If you into it, here's some Mathematica code to check the solution. I believe it's correct but not as close as I would have expected. May be some convergence issues though. nmax = 75; x0 = 0.1; y0 = 0; y1 = 1.; h = 5; m = 4; mysol = NDSolve[{x^2*Derivative[2][y][x] + x*Derivative[1][y][x] + (x^4 + h*x^2 + m^2)*y[x] == 0, y[x0] == y0, Derivative[1][y][x0] == y1}, y, {x, x0, 2}]; sol[x_] := Evaluate[y[x] /. mysol]; p1 = Plot[y[x] /. mysol, {x, x0, 2}]; c1 = I*m; Subscript[a, 0] = 1; Subscript[a, 1] = 0; Subscript[a, 2] = ((-h)*Subscript[a, 0])/((2 + c1)*(1 + c1) + 1 + m^2); Subscript[a, 3] = 0; Table[Subscript[a, n] = -(Subscript[a, n - 4] + h*Subscript[a, n - 2])/((n + c1)*(n + c1 - 1) + 1 + m^2), {n, 4, nmax}]; f1[x_] := Sum[Subscript[a, n]*x^(n + c1), {n, 0, nmax}] f1d[x_] = D[f1[x], x]; c2 = (-I)*m; Subscript[b, 0] = 1; Subscript[b, 1] = 0; Subscript[b, 2] = ((-h)*Subscript[b, 0])/((2 + c2)*(1 + c2) + 1 + m^2); Subscript[b, 3] = 0; Table[Subscript[b, n] = -((Subscript[b, n - 4] + h*Subscript[b, n - 2])/((n + c2)*(n + c2 - 1) + 1 + m^2)), {n, 4, nmax}]; f2[x_] := Sum[Subscript[b, n]*x^(n + c2), {n, 0, nmax}]; f2d[x_] = D[f2[x], x]; myks = First[NSolve[{y0 == k1*f1[x0] + k2*f2[x0], y1 == k1*f1d[x0] + k2*f2d[x0]}, {k1, k2}]]; myy[x_] := k1*f1[x] + k2*f2[x] /. myks; myd1[x_] = D[myy[x], x]; myd2[x_] = D[myy[x], {x, 2}]; N[x^2*myd2[x] + x*myd1[x] + (x^4 + x^2*h + m^2)*myy[x]] /. x -> 0.334 p2 = Plot[myy[x], {x, x0, 2}, PlotStyle -> Red] Show[{p1, p2}]