I have read a paper which drives a recursion relation as follows and starting from Eq.(1) $$ [- \frac{d^2}{dr^2}+\frac{l(l+1)}{r^2}+\frac{\omega^2}{4}r^2+\frac{\zeta}{r}]\,\phi_{nl}(r)=E^r_{nl}\,\phi_{nl} (r) \tag{1}$$
More technically this is the radial Schrodinger equation for a system confined in Harmonic oscillator well. It has been proved that asymptotically correct solutions of the Eq. (1) have the form $$ \phi_nl(r)=r^{l+1}e^{-\omega r^2/4}\, P^{(p)}_{nl}(r) \tag{2} $$ where $\omega>0$ and $$ \lim_{r\to\infty}e^{-\omega r^2/4}\, P^{(p)}_{nl}(r)=\lim_{r\to0} r^{l+1} P^{(p)}_{nl}(r)=0 \tag{3} $$ There exist analytical solutions of Eq.(2) with $$ P^{(p)}_{nl}(r)=\sum^p_{i=0}a_i\,r^i \tag{4} $$ for a discrete set of $\omega_r (p,l)$ and $\zeta_r(p,l)$. By substitution of Eq.(2) and (4) to (1) one gets the recursion relation $$ A_ia_i+B_{i+1}a_{i+1}+C_{i+2}a_{i+2}=0 \tag{5} $$ where $$ A_i=E-\omega_r(i+l+\frac32) \tag{6}\\ B_i=-\zeta_r \\ C_i=i(i+2l+1) $$
Now I want to know how the recursion relation is derived by substitution of Eq (2) and (4) into Eq (1)?
More Precisely: If I add some constants to Eq.(1) as follows how the result (Eqs.(5) and (6)) changes? $$ [- \frac1{2\mu}\frac{d^2}{dr^2}+\frac{l(l+1)}{2\mu r^2}+\frac{\mu\,\omega^2}{4}r^2+\frac{\zeta}{r}]\,\phi_{nl}(r)=E^r_{nl}\,\phi_{nl} (r) \tag{7}$$
When you substitute in (1) you get a second order o.d.e which can be solved by the power series solution method, where you assume that the solution is in the form of a power series polynomial with arbitrary coefficients, and by substitution in (1) you can equate the coefficients of each term of degree say k on the LHS with the corresponds coefficients on the RHS. Then you can find a recursion formula which describe this relationship between the coefficients.