Is there a b̶e̶t̶t̶e̶r̶ generic way to solve linear difference equations with polynomial coefficients?

Question

I was working on a physics problem and came across a linear difference equation of the form $$ a(n,m) f(n+1,m) = b(n,m)f(n,m+2) + c(n,m) f(n, m-2) $$

Where $a,b,c$ are multivariate polynomials in two variables and $f$ is the function to be solved for.

I realized I don't have any good tools / strategies for solving such multivariate linear equations generically. The current strategy that I could come up with is extremely specific to my particular equation in question. My question does anyone know of a generic strategy, similar to series techniques one might use for linear PDEs to write down solutions to the aforementioned class of equations?

for those curious my exact problem is listed below

Solving

$$i\hbar(n+1)a(n+1,m) = - \frac{\hbar^2}{2s}(m+2)(m+1)a(n,m+2) + \frac{1}{2} ka(n,m-2)$$

Where $\hbar, k, s$ can be assumed to constants.

Attempted Strategy:

1. Create a table of functions

$$ \begin{matrix} \vdots & \vdots & \vdots & ⋰ \\ a(0,2) & a(1,2) & a(2,2) & \dots \\ a(0,1) & a(1,1) & a(2,1) & \dots \\ a(0,0) & a(1,0) & a(1,1) & \dots\end{matrix} $$

Observe the table above gives a "linear" relationship between the terms, so that if I fix values for columns $0,1,2,3$ I can then use the equation to autofill terms throughout the entire table

2. Hope to God I can then guess a closed form after staring at the table long enough

Motivation for the Problem (Optional):

When attempting to derive the solution to the quantum harmonic oscillator one writes down the schrodinger equation (we let mass = s as opposed to m, you'll see why later)

$$ E \Psi = H \Psi $$

Recall that total energy is kinetic + potential, i.e. $E = \frac{p^2}{2m} + V$ and furthermore recall that our potential is $\frac{1}{2}kx^2$ in this case where $k$ is the spring constant.

Substituting terms we then have

$$ i \hbar \frac{\partial \psi}{\partial t} = - \frac{\hbar^2}{2s} \frac{\partial^2 \psi}{\partial x^2} + \frac{k}{2}x^2 \psi $$

We can assume $\psi$ is analytic and therefore $\psi(x,t) = \sum_{n,m=0}^{\infty} a(n,m)t^n x^m $

We can now derive the following relationship between the coefficients $a(n,m)$ by substituting the series definition into the schrodinger equation to yield:

$$i\hbar(n+1)a(n+1,m) = - \frac{\hbar^2}{2s}(m+2)(m+1)a(n,m+2) + \frac{1}{2} ka(n,m-2)$$

Answer 1

Well, the canonical procedures would all but take you to the generating function for the a(n,m) s, satisfying the original separable PDE, in so many words, $$\left ( ℏ∂_t +{ℏ^2\over 2}∂^2_x -{\over 2}^2\right )\psi (x,t)=0,$$ as suggested in the comment by @user619894 ; which is easier to solve for: it is obvious that the eigenvalues of $ℏ∂_t$ must cancel versus the eigenvalues of ${ℏ^2\over 2}∂^2_x -{\over 2}^2$.

That is, solutions for have the separated form $\xi_N(t) \phi_N(x)$, with $$ (ℏ∂_t -\lambda_N)\xi_N(t)=0, \qquad \left ({ℏ^2\over 2}∂^2_x -{\over 2}^2 +\lambda_N\right )\phi_N(x)=0, $$ for some spectral index N , not the n of your question; we'll focus on integral N s below.

You then simply observe that the eigenfunctions of the second equation, $$\phi_N=\exp(x^2 \sqrt{ks} /2\hbar)~\partial_x^N ~ \exp(-x^2 \sqrt{ks} /\hbar) $$ for integers N correspond to essentially integrally spaced eigenvalues, $$\lambda_N=\hbar\sqrt{k/s}~(1/2+N) ~.$$

These, then, correspond to exponential eigenfunctions of the first equation, $\xi_N=e^{-i\lambda_N t/\hbar} $.

In all, you have $\psi_N\propto e^{-i \sqrt{k/s}~ N t } ~\exp(x^2 \sqrt{ks} /2\hbar)~\partial_x^N ~ \exp(-x^2 \sqrt{ks} /\hbar) ,$ the generating function of your target matrix elements.