# Solving recursion relation from power series

I am interested in solving differential equations in the form of power series. Let's say we have following equation:

$$f^{\prime \prime} (\rho) + \left( \frac{2 e^{-k \rho}}{\rho} - \varepsilon \right) f(\rho) = 0$$

f''(rho) + (2 Exp(-k*rho)/rho - epsilon) f(rho) = 0   (1)


We can write f(rho) in the form of power series:

f(rho) = Sum[c[n] rho^n, {n,0,max}]
Exp[-k*rho] = Sum[(-k*rho)^n/n!, {n,0,max}]


We can collect terms next to 1/rho, 1, rho, rho^2, ... to get equations for coefficients c, c, ...:

expression = Collect[Sum[c[n + 2] (n + 1) (n + 2) rho^n, {n, 0, max}] +
2/rho Sum[(-k rho)^n/n!, {n, 0, max}] Sum[
c[n] rho^n, {n, 0, max}] - epsilon Sum[
c[n] rho^n, {n, 0, max}], rho]


After specifying c = 0 and c = 1 (c = 0 is required by equations, c = 1 is arbitrary) we get set of equations (from fact that expression = 0, so every coefficient in power expansion have to be zero):

s3 = Solve[(k^2 c - 2 k c - \[Epsilon] c + 2 c + 6 c ==
0) /. {c -> 0, c -> 1, c -> -1}, c];

s4 = Solve[(-(1/3) k^3 c + k^2 c - 2 k c - \[Epsilon] c +
2 c + 12 c == 0) /.
Flatten@{c -> 0, c -> 1, c -> -1, s3}, c];


etc.

Is there a way to automatize this procedure, so that the input would be an equation, or some expression in the form of power series and desired order N and the output would be solution for coefficients c, ..., c[N]?

I thought about something like this:

rho*expression/.rho->0


which gives coefficient next to 1/rho (and kills every other order) but this gives only the lowest order and again I would have to set up the equations manually.

Thanks.

Probably something like this can be done:

n = 5; (* highest-order term *)
f = Sum[C[k] ρ^k, {k, 1, n}] + O[ρ]^(n + 1); (* C == 0 already imposed *)

Solve[Thread[CoefficientList[D[f, {ρ, 2}] + (2 Exp[-k ρ]/ρ - ε) f, ρ] == 0],
Array[C, n]] // Simplify
{{C -> -C, C -> 1/6 (2 + 2 k + ε) C,
C -> -(1/36) (2 + 8 k + 3 k^2 + 4 ε) C,
C -> 1/360 (2 + 33 k^2 + 6 k^3 + 10 ε + 3 ε^2 + 4 k (5 + 3 ε)) C}}


along with a Solve::svars warning since C was not specified. Increase n as seen fit.

eq= f''[x]+(2 Exp[-k*x]/x-e) f[x]==0


For a given value of max, you are looking for a series solution sol:

max=5; sol=Function[x, Sum[c[i] x^i, {i,1,max}] + O[x]^max]


Substitute this solution in your equation and use the function SolveAlways for finding the coefficients:

SolveAlways[eq /. f->sol, x]

(* {{c->-(1/36) (2+4 e+8 k+3 k^2) c,c->1/6 (2+e+2 k) c,c->-c}} *)

• Now I don't know which solution should I accept, both of them do the job, but this one is perhaps a bit more clear... If I try this for max at least 7 I get two sets of solutions, one with "correct" (non-zero) coefficients, the second with all coefficients being zero. Why is that so? – user16320 Apr 11 '16 at 22:33
• It seems that for max at least 7, the function SolveAlways does not give correct results, A rule for the variable e is returned as well. Therfore, I would accept the other, correct, answer. – Fred Simons Apr 12 '16 at 7:34