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.