Sage code:
values = {}
def A(n,k):
if (n,k) in values:
return values[(n,k)]
if n==0:
if k==0:
result = 1
else:
result = 0
else:
result = (n + 1 - k)**2 * A(n-1, k-1)
result += 2*(n*(n+1)-k**2) * A(n-1, k)
result += (n + 1 + k)**2 * A(n-1, k+1)
values[(n,k)]=result
return result
for n in range(5):
print [A(n,k) for k in range(-n, n+1)]
Output:
[1]
[1, 4, 1]
[1, 26, 66, 26, 1]
[1, 120, 1191, 2416, 1191, 120, 1]
[1, 502, 14608, 88234, 156190, 88234, 14608, 502, 1]
These are famous numbers to combinatorists. So it looks like you're getting Eulerian numbers for odd n, with k shifted. See A008292. Note that it's easy to see that A(n,k) = A(n,-k) by induction on n, even if you didn't know that.