Let $S_k^n$ be the number of possible surjections from a set of $k$ elements to a set of $n$ elements. We have $$\begin{align} &S_0^0 = 1,\qquad\forall k>0: S_k^0 = 0,\\ &S_n^n = n!,\qquad\forall k<n: S_k^n = 0,\\ &\forall k>n>0: S_k^n = n(S_{k-1}^n + S_{k-1}^{n-1}). \end{align}$$ The last relation can be seen as follows: let $\phi:\underline{k}\to\underline{n}$ be a surjection. Consider the subset $\underline{k-1}\subset\underline{k}$ consisting of the first $k-1$ elements $\{1,2,\ldots,k-1\}$. Then either the restriction of $\phi$ to $\underline{k-1}$ is a surjection, and $\phi(k)$ is any element of $\underline{n}$, or the restriction misses exactly one element among the $n$ elements of $\underline{n}$, and $k$ is mapped to this missed element. Representing the first terms in a grid, we have $$\begin{array}{c|cccccc} n\backslash k & 0 & 1 & 2 & 3 & 4 & 5 & 6\\ \hline 0 & 1 & \color{lightgray}{0} & \color{lightgray}{0} & \color{lightgray}{0} & \color{lightgray}{0} & \color{lightgray}{0} & \color{lightgray}{0} \\ 1 & \color{lightgray}{0} & 1 & 1 & 1 & 1 & 1 & 1\\ 2 & \color{lightgray}{0} & \color{lightgray}{0} & 2 & 6 & 14 & 30 & 62\\ 3 & \color{lightgray}{0} & \color{lightgray}{0} & \color{lightgray}{0} & 6 & 36 & 150 & 540\\ 4 & \color{lightgray}{0} & \color{lightgray}{0} & \color{lightgray}{0} & \color{lightgray}{0} & 24 & 240 & 1560 \end{array}$$ If we extend $S_k^n$ to $n,k<0$ by setting it to zero everywhere, then the last relation is satisfied everywhere except at the point $(k,n) = (0,0)$.
I am interested in the generating function of these numbers, that is $$S(x,y) = \sum_{k,n\in\mathbb{Z}} S_k^nx^ky^n.$$ Using the relation above, we obtain $$\begin{align} S(x,y) = & 1 + \sum_{k,n\in\mathbb{Z}}n(S_{k-1}^n + S_{k-1}^{n-1})x^ky^n\\ = & 1 + xy\left(S(x,y) + (1+y)\frac{\partial S}{\partial y}(x,y)\right). \end{align}$$ Moreover, we have the boundary conditions $S(x,0) = 1 = S(0,y)$. According to WolphramAlpha, this is solved by $$S(x,y) = c(x)\left(\frac{y}{(y+1)^{x+1}}\right)^{\frac{1}{x}} + (y+1)^{-\frac{x+1}{x}}{}_2F_1\left(-\frac{1}{x},-\frac{1}{x};\frac{x-1}{x};-y\right).$$ (This was edited after an error in the original question was found by @MickA)
Now, just by looking at the equation we know that this function will automatically satisfy the required boundary conditions. How can we fix the function $c(x)$? I am a bit at loss (also because I'm not familiar with hypergeometric functions at all). Also, is it possible derive a closed expression for $S_k^n$ from the resulting solution?
Remark: I already know of the formula obtained by inclusion-exclusion principle: $$S_k^n = \sum_{i = 0}^n(-1)^i\binom{n}{i}(n-i)^k$$ for $0\le n\le k$. I am curious to see if this alternative method works, if it gives exactly the same answer, and if it contains any additional information in general.