I have a two-variable recurrence relation of the form,
\begin{align} -&[(N+1)n+N(n+1)+(M+1)m+M(m+1)]p(n,m)\\ -&\epsilon[(n+1)m+(m+1)n)]p(n,m)\\ +&(N+1)(n+1)p(n+1,m)+(M+1)(m+1)p(n,m+1)\\ +&N n p(n-1,m)+M m p(n,m-1)\\ +&\epsilon n(m+1)p(n+1,m-1)+\epsilon m(n+1)p(n-1,m+1)\\ =&0, \end{align}
where $p(n,m)$ are probabilities, i.e., $p(n,m)\geq0 ~ \forall n,m $ and $\sum_{n=0}^\infty\sum_{m=0}^\infty p(n,m)=1$. I also emphasize that $n,m$ take values from $0$ to $\infty$, and then $p(-n,m)=0=p(n,-m)$ is required.
I've been struggling for some time in solving this problem, which stems from a physical model. In this regard, I must say that I have already numeric solutions and am now looking for an analytic solution. Yet, in my physical model $\epsilon<<1$, and a solution at first order in $\epsilon$ is welcome.
Moreover, I have the solution for $\epsilon=0$; in this regime $n$ and $m$ are independent and \begin{align} p_0(n,m)= C \left(\frac{N}{N+1}\right)^n\left(\frac{M}{M+1}\right)^m, \end{align} with a constant $C$ fixed by normalization of the probabilities.
Since I have little experience with recurrence relations, my tentative solution is to map this problem into differential equations.
My trial
I consider the Z-transform of the probability $p(n,m)$,
\begin{align} F(u,v):=\mathcal{Z}[p]=\sum_{n,m=0}^{\infty} u^{-n} v^{-m} p(n,m), \end{align} whose inverse is given by the contour integral in the two-dimensional complex-plane \begin{align} p(n,m)=\frac{1}{(2\pi i)^2 }\oint_{C_1}\oint_{C_2} u^{n-1} v^{m-1} F(u,v), \end{align} with $C_1, C_2 $ counterclockwise and enclosing all singularities of the integrated.
Applying the transform to the recurrence relation, I have the following PDE for $F(u,v)$,
\begin{align} (\mathcal{D}_0+ \epsilon \mathcal{D}_1)F=0, \end{align} with \begin{align} \mathcal{D}_0&= (u^{-1}-1)(N+ u(u(N+1)-N)\partial_u)\\ +&(v^{-1}-1)(M+ v(v(M+1)-M)\partial_v)\\ \mathcal{D}_1&=(u-v)^2 \partial_{uv} + u(1-v^{-1})\partial_u + v(1-u^{-1})\partial_v. \end{align}
As before, I have the $\epsilon=0$ solution, \begin{align} F_0(u,v)= \frac{u}{u(N+1)-N}\frac{v}{v(M+1)-M}, \end{align} satisfying $\mathcal{D}_0 F_0 =0$. To enforce normalization, I have imposed that $F_0(1,1)=1$.
The full PDE seems to be complicated, and I am not really sure about the correct boundary conditions. My guess would be that they should come from requiring that \begin{align} p(-1,m)=\frac{1}{(2\pi i)^2}\oint_{C_1}\oint_{C_2} u^{-2} v^{m-1} F(u,v)=0,\\ p(n,-1)=\frac{1}{(2\pi i)^2}\oint_{C_1}\oint_{C_2} u^{n-1} v^{-2} F(u,v)=0, \end{align} but I am not really sure. I think if I get this boundary conditions right, Mathematica can solve the problem at first order in $\epsilon$. This would boil down to assuming that $F=F_0+\epsilon \phi$, which implies
\begin{align} \mathcal{D}_0\phi = - \mathcal{D}_1F_0 + \mathcal{O}(\epsilon^2). \end{align} That is, the problem is to solve the inhomogeneous equation for $\phi$, where $\phi(1,1)=0$, since the inverse-transform of $\phi$ must preserve the normalization of $F_0$, i.e. $F(1,1,)=1$.
Maybe this is not the easiest way to solve the problem but seemed the most familiar to me. Any help on my trial or in the original statement of the problem is very welcome. To do it pen and paper I think the path should be to find the Green's function for $\mathcal{D}_0$, but the lack of boundary conditions is still a problem. Thanks in advance.