Finding a solution to this equation using Frobenius method

[climbon]

Hi, I have this equation to solve.

y'' + (1/x)y' + [(x^2) + k + (m^2 / x^2)]y = 0

now, I've tried to solve this using frobenius method but cannot formulate a solution.

I have that

a_(n+4) = [-ka_(n+2) - a_(n)] / [n^2 +/- 2inm]

is my recurrence relation, but now I'm stuck and don't know how to proceed, any help greatly appreciated.

thanks.

 

[jackmell]

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}]