Let $\beta >0 $. The question is to solve a following recursion: \begin{equation} P^{(j+2)}(\beta) = \frac{\imath}{2} \left[ \left((-1+\beta) j - 1\right) P^{(j+1)}(\beta) + \sum\limits_{j_1=1}^{j+1} P^{(j_1)}(\beta) P^{(j+2-j_1)}(\beta)\right] \end{equation} for $j=-1,0,1,2,\dots$. The initial conditions are $P^{(0)}(\beta) = \imath$ and $P^{(1)}(\beta) = \beta/2$.
It is customary to use the generating functions technique to solve the recursion. We define a generating function as: \begin{equation} F(x) := \sum\limits_{j=0}^\infty P^{(j)}(\beta) x^j \end{equation}
Now we carry out a Z-transform. We multiply both sides of our original equation by $x^{j+2}$ and and sum over $j=-1,0,1,\dots$. We have:
\begin{equation} -P^{(0)} + F(x) = \frac{\imath}{2}\left( (-1+\beta) \sum\limits_{j=-1}^\infty j P^{(j+1)} x^{j+2} - x F(x) \right) + \frac{\imath}{2} \left(F(x) - P^{(0)}\right)^2 \end{equation} The sum in the brackets on the right-hand-side is equal to $x\left(x F^{'}(x) - F(x)\right)$. After simplifying the whole expression we get a first order differential equation for the generating function . We have:
\begin{equation} 1 - \beta x F(x) + F(x)^2 = (1-\beta) x^2 F^{'}(x) \end{equation}
Before we proceed we do a quick sanity check. Setting $x=0$ we get $F(0) = \pm \imath$ which matches the initial condition $P^{(0)} = \imath$. Now, back to the differential equation. We recognise that it is a Riccati equation and Wikipedia tells us that as such it can be reduced to a linear second order differential equation. We use Mathematica to find the solution. We have: \begin{equation} F(x) = \frac{J_{-1 + \frac{1}{2(1-\beta)}}\left[\frac{1}{(-1+\beta) x}\right] \Gamma\left[\frac{1}{2(1-\beta)}\right] + C J_{1 - \frac{1}{2(1-\beta)}}\left[\frac{1}{(-1+\beta) x}\right] \Gamma\left[\frac{1}{2(-1+\beta)}\right]} {-J_{\frac{1}{2(1-\beta)}}\left[\frac{1}{(-1+\beta) x}\right] \Gamma\left[\frac{1}{2(1-\beta)}\right] + C J_{\frac{1}{2(-1+\beta)}}\left[\frac{1}{(-1+\beta) x}\right] \Gamma\left[\frac{1}{2(-1+\beta)}\right]} \end{equation}
Now, we plot the function above for $\beta= 0.9000001$ and $C = 0.5$.
As we can see that function has infinitely many singularities in the vicinity of the origin. This means that we cannot even compute $F(0)$ not to mention the higher derivatives.On the other hand we need to compute those derivatives since the inverse of the Z-transform is given by:
\begin{equation} P^{(j)}(\beta) = \left. \frac{1}{j!} \frac{d^j F(x)}{d x^j} \right|_{x=0} \end{equation}
However, since our generating function is singular at the origin there is no way for us to to invert the Z-transform. Can anybody tell me where is a bug in this whole reasoning? How do we invert the Z-transform ?
Now, let me mention where did I come across this problem. Originaly I was trying to solve a following differential equation: \begin{equation} (I) \frac{1}{\omega^2} \frac{d^2 r_\xi}{d \xi^2} + \xi^{-2 \beta} r_\xi = 0 \end{equation} using the Wentzel–Kramers–Brillouin (WKB) approximation. This approximation holds when $\omega$ is ``big''. In this approximation one assumes that the solution is given as: \begin{equation} (II) r_\xi = \exp\left[\omega \sum\limits_{j=0}^\infty \frac{S_j(\xi)}{\omega^j} \right] \end{equation} Insering the ansatz (II) into the equation (I) and equating the terms at consecutive powers of $\omega$ to zero we get a series of equations for the functions $S_j(\xi)$. Let is relatively easy to write a Mathematica program that performs this procedure. Having done that we discover that the functions $S_j(\xi)$ are given by:
\begin{equation} \dot{S}_j(\xi) = P^{(j)}(\beta) \xi^{(-1+\beta) j - \beta} \end{equation} for $j=1,2,\dots$ and $\dot{S}_0(\xi) = \imath \xi^{-\beta}$. It appears that the polynomials $P^{(j)}(\beta)$ satisfy the recursion relations in the formualtion of the problem.