We study the concentration of a degree-$d$ polynomial of the $N$ spins of a general Ising model, in the regime where single-site Glauber dynamics is contracting. For $d=1$, Gaussian concentration was shown by Marton (1996) and Samson (2000) as a special case of concentration for convex Lipschitz functions, and extended to a variety of related settings by e.g., Chazottes et al. (2007) and Kontorovich and Ramanan (2008). For $d=2$, exponential concentration was shown by Marton (2003) on lattices. We treat a general fixed degree $d$ with $O(1)$ coefficients, and show that the polynomial has variance $O(N^d)$ and, after rescaling it by $N^{-d/2}$, its tail probabilities decay as $\exp(- c\, r^{2/d})$ for deviations of $r \geq C \log N$.