Solving recursion relation from power series

Question

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.

Answer 1

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.

Answer 2

This is your equation:

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