Closed form of a sequence defined by two indices and a recursion formula Let $N$ be a strictly positive integer.
Define the following sequence $u$ for integers $p,q \in [1,N] \times \mathbb{Z}$:
$u(p,q) = 0$ if $q<0$ or $q>N$.
$u(1,q) = \delta(1,q)$ ( $1$ if $q=1$ else $0$)
$u(p,q) = \frac{N-q+1}{N}u(p-1,q-1)+\frac{q+1}{N}u(p-1,q+1)$
Question:
I would like to find a closed form for $u(N,0)$.
We have $u(N,0) = \frac{1}{N}u(N-1,1)=\frac{N-1}{N}u(N-2,0)+\frac{2}{N}u(N-2,2) = ...$
We can apply iteratively the recursion formula until we attain the initial conditions.
But I don't know if it is possible to obtain a closed-form for this?
Obviously if we can find a closed-form for $u(p,q)$ then we are done.
Any hints are appreciated. Thank you.
 A: I don't know if there is a non-recursive formula for $u(N,0)$, and in this answer we only reduce the original problem to the problem of finding a (non-recursive) solution of a system of equations given by a (relatively simple) tridiagonal matrix.
There is an underlying one-dimensional Markov system. Consider the states $\{0,\dots,N\}$, and the transition probabilities given by $P(i\to i+1)= \frac{N-i}{N}$, $P(i\to i-1)=\frac iN$ and $P(i\to j)=0$ if $|i-j|\ne 1$. Then the transition matrix is the tridiagonal matrix
$$
M=\begin{pmatrix}
0&\frac 1N &0&\dots&\dots&0&0\\
1&0&\frac 2N&\dots &\dots &\dots &0\\
0&\frac{N-1}{N}&0&\frac 3N&&&\vdots\\
\vdots&&\ddots&&\ddots &&\vdots\\
\vdots&&&\ddots&&\frac{N-1}{N}&0\\
\vdots &&&&\ddots &0&1\\
0&\dots&&&&\frac 1N&0
\end{pmatrix}
$$
Then $u(i,j)=(M^i)_{j,0}$, in particular $u(N,0)=(M^N)_{0,0}$. 
Although it is written in a short formula, one may still consider this as a recursive formula.
If we could know the decomposition of $M$ as $PDP^{-1}$, where $D$ is the diagonal matrix with the eigenvalues, then we would obtain a closed formula for $u(N,0)$, since
then 
$$
u(N,0)=(M^N)_{0,0}=(PD^NP^{-1})_{0,0}.
$$
With the help of Mathematica we obtain the following results (without proof):


*

*The matrix $D$ is known when $N$ is even. In that case the $N$ eigenvalues of $M$ are $\left\{\pm \frac{2k}N\right\}_{k=1,\dots,N/2}$.

*We also know that the eigenvector corresponding to $\lambda=1$ is $v=(v_i)$ with $v_i=\binom{N}{i}$. 

*If $w=(w_0,w_1,\dots,w_N)$ is an eigenvector corresponding to the value $\lambda=\frac{2k}{N}$, then the vector given by $\bar w=((-1)^iw_i)_i$ is an eigenvector corresponding to
$\lambda=-\frac{2k}{N}$. 
Here is another approach:
Usually one tries to write $u(k,0)$ as $A_1\lambda_1^k+A_2\lambda_2^k+\dots A_N\lambda_N^k$
for some coefficients $A_i$ and the (known) eigenvalues $\lambda_i$, but in order to obtain 
a system of $N$ equations in order to determine the $A_i$'s, one need $N$ steps, so it wouldn't help if you want to know $u(N,0)$. 
On the other hand the asymptotics may provide some condition on $A_1,A_2$, where $\lambda_1=1$ and $\lambda_2=-1$, so maybe we only need $N-2$ steps.
Since the columns of the matrix $P$ are the eigenvectors, the problem is reduced to the following:
${\bf Problem:}$
For $k=1,\dots,N/2$ find a (non-recursive) solution $v$ to the system
$$
(1)\quad\quad\quad \quad\quad\quad \quad\quad\quad Mv=\pm\frac{2k}N v.
$$
This is equivalent to determining the Kernel of the tridiagonal matrix
$$
(M\mp\frac{2k}N Id)=\begin{pmatrix}
\mp\frac{2k}N&\frac 1N &0&\dots&\dots&0&0\\
1&\mp\frac{2k}N&\frac 2N&\dots &\dots &\dots &0\\
0&\frac{N-1}{N}&\mp\frac{2k}N&\frac 3N&&&\vdots\\
\vdots&&\ddots&&\ddots &&\vdots\\
\vdots&&&\ddots&&\frac{N-1}{N}&0\\
\vdots &&&&\ddots &\mp\frac{2k}N&1\\
0&\dots&&&&\frac 1N&\mp\frac{2k}N
\end{pmatrix}.
$$
I don't know if there is a (non-recursive) solution $v$ for the system (1).
${\bf Edit:}$
I found $A_i=\frac{1}{2^N}\binom{N}{N/2-k}$ for $\lambda_i=\pm \frac{2k}{N}$. Hence, if $N=2m$, we have
$$
(2)\quad\quad \quad\quad \quad\quad \quad\quad \quad\quad u(N,0)=\frac{2}{N^N}\sum_{k=1}^{m}\binom{N}{m-k}k^{N}=\frac 1{N^N}\sum_{k=0}^{N}(m-k)^{N}\binom{N}{k}.\quad\quad \quad\quad \quad\quad \quad\quad 
$$
${\bf Second\ Edit\ (Short\ Answer):}$
There is an interpretation of $N^N u(N,0)$ as the number of functions $f:\{1,\dots,N\}\to \{1,\dots,N\}$ such that every preimage has an even cardinality.
From this interpretation it follows that $N^N u(N,0)=a(N/2)$, where $a(n)=$https://oeis.org/A209289
For  $N$ odd we have $u(N,0)=0$. So assume $N=2m$. Write $\Bbb{N}_i=\{1,\dots,i\}$. Define $U(i,j)$ to be the number of functions from $\Bbb{N}_i$ to $\Bbb{N}_N$ such that there are $j$ elements in $\Bbb{N}_N$ with a preimage that has an odd cardinality, i.e.,
$$
U(i,j)=\#\Big\{f:\Bbb{N}_i\to \Bbb{N}_N\ \big|\ j=\#\{x\in\Bbb{N}_N,\ | \ f^{-1}(x) \ odd\}\Big\}.
$$
Then 
$$
U(i,j)=(j+1)U(i-1,j+1)+(N-j+1)U(i-1,j-1).
$$
Using this recursive relation it is easy to see that $\frac{1}{N^i}U(i,j)=u(i,j)$. In  particular 
$$
u(N,0)=\frac 1{N^N}U(N,0).
$$
So $N^N u(N,0)$ is the number of functions $f:\{1,\dots,N\}\to \{1,\dots,N\}$ such that every preimage has an even cardinality, which is $a(N/2)$ for $a(m)$ given by 
https://oeis.org/A209289, and so $u(N,0)$ is given by (2).
