For $k=1,\ldots, \infty$, and $n=2,\ldots, k-1$ we define the following recurrence relation $$ p_{1,k} = \frac{k-1}{k} p_{1,k-1}$$ $$ p_{n,k} = \frac{k-1}{k} p_{n,k-1}+\frac{1}{k^2} p_{n-1,k-1}$$ $$ p_{k,k} = \frac{1}{k^2} p_{k-1,k-1}$$ $$ p_{1,1} = 1$$
We may evaluate $p_{1,k}$ to be given by $$ p_{1,k} = \frac{1}{k}$$
Similarly, we may evaluate $p_{2,k}$ to be given by $$ p_{2,k} = \frac{k-1}{k^2}$$
We seek a generalised relation for $p_{n,k}$ which is defined non-recursively. Is this possible? If so, what form may this take?
A common approach to these kinds of problem is through the use of generating functions.
I'll show the start of one possible approach to find these coefficients using a generating function, but I haven't worked all the way through to determine whether it's the best approach.
First, we'll note that $p_{k, k} = \frac{1}{(k!)^2$, which we'll make use of. So overall our recursion has become:
$$\begin{eqnarray} p_{1, k} & = & 1 \\ p_{k, k} & = & \frac{1}{(k!)^2} \\ p_{j, k} & = & \frac{1}{k^2} p_{j-1, k-1} + \frac{k - 1}{k} p_{j, k-1}\ &\mbox{for } 2 \leq j \leq k - 1 \end{eqnarray}$$
So now we can define our generating function. It will be defined as:
$$A(x, y) = \sum_{k = 1}^\infty \sum_{j = 1}^k p_{j, k} x^j y^k$$
We can then expand that out and apply our recursion to get something like:
$$\begin{eqnarray} A(x, y) & = & \sum_{k = 1}^\infty \sum_{j = 1}^k p_{j, k} x^j y^k \\ & = & x y + x y^2 + \frac{1}{4} x^2 y^2 + \sum_{k = 3}^\infty \left[x y^k + \frac{1}{(k!)^2} (xy)^k + \sum_{j = 2}^{k - 1} \left(\frac{1}{k^2} p_{j - 1, k - 1} + \frac{k - 1}{k} p_{j, k - 1} \right) x^j y^k \right] \end{eqnarray}$$
You can then break up the last term into a number of separate terms, and through some careful rearranging get some terms involving $p_{j, k} x^j y^k$ that you can hopefully transform back into a version of the original function (noting that if you see terms like $k p_{j, k} y^k$ then that's indicative of derivatives).
Some other potentially useful identities are $\sum_{k = 1}^\infty z = \frac{z}{z - 1}$ and $\sum_{k = 0}^\infty \frac{z^k}{(k!)^2} = I_0 (2 \sqrt{z})$, where $I_0(z)$ is a modified Bessel function of the first kind.
Given the structure of the recurrence, making the generating function exponential in $y$ might be more useful, or there might be a transformation of the original recurrence that would simplify things (I played around with definitions like $q_{j, k} = \frac{k!}{(k - j)!} p_{j, k}$ and they showed a little promise).