# 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[0], c[1], ...:

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] = 0 and c[1] = 1 (c[0] = 0 is required by equations, c[1] = 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[0] - 2 k c[1] - \[Epsilon] c[1] + 2 c[2] + 6 c[3] ==
0) /. {c[0] -> 0, c[1] -> 1, c[2] -> -1}, c[3]];

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

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[0], ..., 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] == 0 already imposed *)

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

along with a Solve::svars warning since C[1] 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[4]->-(1/36) (2+4 e+8 k+3 k^2) c[1],c[3]->1/6 (2+e+2 k) c[1],c[2]->-c[1]}} *)
• 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? Apr 11, 2016 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. Apr 12, 2016 at 7:34