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.