For a finite non-empty set $\Omega=\{\omega_1,\dots,\omega_{N+1}\}\subset \mathbb{R}$, denote by $\Phi_{\Omega}$ the linear functional which maps a function $g:\Omega\mapsto \mathbb{R}$ to $$\Phi_{\Omega}(g)=\sum_{\omega\in \Omega} \frac{g(\omega)}{\prod_{\tau\in \Omega\setminus \omega} (\omega-\tau)}.$$ This expression may be viewed as a coefficient of $x^{N}$ in the polynomial $P_{\Omega}(g)$ which interpolates $g$ on the set $\Omega$. Note that for any polynomial $p(x)\in \mathbb{R}[x]$ we have $\Phi_{\Omega}(p)=[x^N]P_{\Omega}(p)$ (as usual, $[x^N]$ stands for a coefficient of $x^N$) and $P_{\Omega}(p)$ is a remainder of $p(x)$ modulo the polynomial $\prod_{\omega\in \Omega}(x-\omega)$. Therefore $\Phi_{\Omega}(h)=0$ whenever $\deg h<|\Omega|-1$, and $\Phi_{\Omega}(x^{N+k})=h_k(\Omega)$, where $h_k(\Omega)$ stands for $$h_k(\Omega)=h_k(\omega_1,\dots,\omega_{N+1})=\sum_{i_1\leqslant i_2\dots\leqslant i_k} \omega_{i_1} \omega_{i_2}\dots \omega_{i_k},$$ the complete symmetric polynomial. This is a basic fact from the theory of symmetric functions, the short way to see it is to observe that $h_k(\Omega)$ is a coefficient of $x^{k}$ in the product $\prod_{\omega\in \Omega}(1+\omega x+\omega^2x^2+\dots)=\prod_{\omega\in \Omega} (1-\omega x)^{-1}$ and expanding this rational function as a linear combination of elementary fractions $1/(1-\omega x)$.
You are interested in the case when $\Omega=\{0,1,\dots,N\}$ and $g(x)=x^Nf(x)$. It is less or more natural to expect that $\Phi_{\Omega}(g)$ behaves like the coefficient of $[x^N]$ in the Taylor polynomial $p_N(x)=\frac1{N!}g^{(N)}(N/2)(x-N/2)^N+\dots$ of the function $g$ expanded in the point $N/2$. Let us justify this for $g(x)=x^{N-a}$, this gives your asymptotics $\frac1{N!}g^{(N)}(N/2)=(N!)^{-1}(N/2)^{-a}(1-a)(2-a)\dots(N-a)\sim 2^aN^{-2a}/\Gamma(1-a)$ for large $N$.
We need the following
Lemma. Assume that $\Omega$ is a symmetric w.r.t. 0 set. Then $h_k(\Omega)=0$ for odd $k$ and $h_{2m}(\Omega)=h_m(\Omega^2\setminus 0)$, where $\Omega^2=\{\omega^2:\omega\in \Omega\}$. In particular, this gives the estimate $|h_{2m}(\Omega)|\leqslant \max(\Omega)^{2m}\cdot \binom{[N/2]+m-1}m$.
Proof. The formula is seen from the generating function $\prod 1/(1-\omega x)$. The estimate is obtained by counting the number of terms and estimating each of them as $\max(\Omega)^{2m}$.
Now write $g(x):=x^Nf(x)=p_N(x)+r_N(x)$, we have $\Phi_{\Omega}(g)=\Phi_{\Omega}(p_N)+\phi_{\Omega}(r_N)$, where $p_N(x)=\frac1{N!}g^{(N)}(N/2)(x-N/2)^N+\dots$ is the Taylor polynomial of the function $g(x)$ at the point $N/2$ and $r_N(x)$ is the remainder. As said beforуб the first summand $P:=\Phi_{\Omega}(p_N)$ behaves as expected, and it remains to estimate the second summand. It is more convenient to do it for any $\Omega\subset (0,N)$ symmetric w.r.t. $N/2$, the estimate would be uniform in $\Omega$, so it continues to $\Omega\subset [0,N]$.
You may expand $r_N(x)$ as a sum $\sum_{k\geqslant 1} \frac1{(N+k)!}g^{(N+k)}(N/2)(x-N/2)^{N+k}$. Let $\Omega=\{N/2+\rho_1,N/2+\rho_2,\dots,N/2+\rho_N\}$, we have $\rho_i\in (-N/2,N/2)$ and the set $\Omega_s=\Omega-N/2$ is symmetric w.r.t. 0. Now the series in $k$ and the interpolating operator $\Phi_{\Omega}$ may be permuted (simply because the sum in this operator is finite for any fixed $N$). We get $$\Phi_{\Omega}(r_N)=\sum_{k\geqslant 1} \frac1{(N+k)!}g^{(N+k)}(N/2) h_k(\Omega_s).$$ Comparing $\frac1{(N+k)!}g^{(N+k)}(N/2)$ and $P=\frac1{N!}g^{(N)}(N/2)$ we see (remember that $g^{(N)}(x)=Cx^{-a}$ for certain $C$) that $$ \left|\frac{P^{-1}}{(N+k)!}g^{(N+k)}(N/2)\right|\leqslant \frac{(2/N)^k k! k^{A}}{(N+1)(N+2)\dots (N+k)} $$ where a constant $A$ depends on $a$ (so that $|a|(|a|+1)\dots (|a|+k)\leqslant k! k^A$ for all $k\geqslant 2$).
Multiplying this by the estimate provided by Lemma we totally get (denote $k=2m$) at most $$ \sum_{m=1}^\infty \frac{(2/N)^{2m} (2m)! (2m)^{A}}{(N+1)(N+2)\dots (N+2m)}\cdot \frac{(N/2+1)\dots (N/2+m)}{m!}\cdot (N/2)^{2m} $$ It is easy that if we denote $m$-th term by $x_m$, than $x_m=O_a(m^{-2}N^{-1})$, where a constant in $O(\cdot)$ depends only on $A$ (which in turn depends only on $a$): look at the ratio $x_{m+1}(m+1)^2/(x_m m^2)$, it eventually decreases. This gives $O(1/N)$ after summation by $m$.
(to be continued for other functions $f$).
For $\psi'$ you may use the formula $\psi'(x)x^N=\int_0^\infty \frac{t}{1-e^{-t}}e^{-xt}x^Ndt$ and interchange the integral and interpolating functional, since interpolating the function $x^Ne^{-xt}$ in the points $0,1,\dots,N$ is sort of pleasant: $$ \frac1{N!}\sum_{n=0}^N(-1)^{N-n}\binom{N}ne^{-nt}n^{N}=\frac1{N!}(d/dt)^N (1-e^{-t})^N $$ We integrate by parts to get the answer in the form $$ \int_0^\infty -dt\cdot \frac{d}{dt}\frac{t}{1-e^{-t}}\cdot \frac1{N!}(d/dt)^{N-1} (1-e^{-t})^N. $$ Next, we write $$ \frac1{N!}(d/dt)^{N-1} (1-e^{-t})^N=\frac1N[z^{N-1}](1-e^{-t-z})^N $$ for using Lagrange inversion formula in the following form: $$ \frac1N[z^{N-1}]\varphi(z)^N=[w^N] g(w) $$ where $f(z)=z/\varphi(z)$ [ah, of course, it is not at all your $f$ from the post] and $g$ is functional inverse of $f$. In our situation $f(z)=f_t(z)=z/(1-e^{-t-z})$ and so $g(w)=g_t(w)$ is the solution of the equation $g/(1-e^{-t-g})=w$. The growth of coefficients of the analytic function $g$ is determined by its first singularity. The first singularity of the inverse function $g$ is the first critical value of the initial function $f$ (that is, in the point closest to the origin where $f'=0$). Consider $t=0$, I guess for $t>0$ the singularities are further from 0 (need to check of course). Then $f(z)=z/(1-e^{-z})$, the equation $f'(z)=0$ rewrites as $e^{z_0}=z_0+1$, then $f(z_0)=z_0+1$ is a solution of $x_0=e^{x_0-1}$ as you conjecture.
(to be continued)
