Dnk[n_, k_] := Dnk[n, k] = Which[k < 0, 0, k == 0, Binomial[2*n, n], True, Dnk[n + 1, k - 1] - Dnk[n, k - 1]];
T[n_, k_] := Module[{xyvec, i, x , ip}, xyvec = {}; For[i = 0, i <= n, i++, AppendTo[xyvec, {i, Dnk[i, n]*Product[i + j, {j, 1, n}]/Dnk[i, 0]}]]; ip = InterpolatingPolynomial[xyvec, x]; Coefficient[ip, x, n - k]];