Edit: ok, now that I have more than 5 minutes to spare I can clean this up a bit and add a wikipedia reference.

I'm going to write A(n,k) for $A_n(k)$. First of all, note that it's easy to see that A(n,k) = A(n,-k) by induction on n, and that the A(n,k) are zero unless -n <= k <= n. So we may as well just start computing these things (with dynamic programming, for good practice) before we start thinking terribly hard:

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]
```

One glance at the third row will tell any combinatorist that these are Eulerian numbers (at least, for odd n). See sequence A008292 at oeis.org. Also, wikipedia has a perfectly reasonable page on the Eulerian numbers: http://en.wikipedia.org/wiki/Eulerian_number. There you can find a recursive formula. I'll use E(n,m) since A is taken already:

$E(n,m) = (n-m)E(n-1,m-1) + (m+1)E(n-1,m)$.

Of course this notation is different than yours; I think your numbers are $E(2n+1, m-n)$, You should be able to see this by applying the above recursive formula twice and doing the above change of variables to recover your own formula, though I haven't done it and may have made an error. There's lots of formulas for the Eulerian numbers and there's a lot known about them.