I'm interested in solving a particular non-linear recurrence in two variables:

$$\lambda_{j,k} = (j+k) \lambda_{j,k-1} + \begin{pmatrix} j+k \\ 2 \end{pmatrix} \lambda_{j-1, k-1}.$$

Here $j \geq -1$ and $k \geq 0$, and we have initial conditions:

$\lambda_{-1,k} = 0$ for all $k$;

$\lambda_{j,0} = 0$ for all $j>0$;

$\lambda_{0,0} = 1$.

This is a relation between the leading coefficients of certain polynomials occurring in a problem in modular invariant theory. $\lambda_{-1,k}$ doesn't really make sense in context, but introducing it simplifies the relation and initial conditions nicely. Because of the context, I happen to know that for any prime $p$, $\lambda_{j,k} \neq 0$ mod $p$ unless $j+k \geq p$. This suggests to me that there is a particularly simple solution, possibly a product of binomial coefficients.

Some starting points: Obviously $\lambda_{0,k} = k!$ and it's easy to show that $\lambda_{j,k} = 0$ if $j>k$ and $\lambda_{k,k} = \frac{(2k)!}{2^k}$ for $k \geq 0$. I've tried solving it with generating functions in two variables, but that just produces a truly terrifying PDE of order 2 and degree 4.