Why should we suspect that the recurrence $T(n) = T(n-1) + n(n-1)$ satisfies a polynomial identity? In the question Algorithms: Recurence Relation, the author asked about the recurrence relation $$T(n) = T(n-1) + n(n-1)$$ and one of the answers proposed assuming $T(n)$ is polynomial, then manipulating the equation to obtain the coefficients.
Question: What reason is there to suspect that this recurrence satisfies a polynomial identity?  And how can we apply this to linear recurrences in general?
The Fibonacci Numbers, for example, instead satisfy a exponential solution.  So, if we attempted to fit a polynomial solution, it wouldn't work.
 A: Let $a_0,a_1,a_2,a_3,\dots$ be a sequence.  Let $b_0,b_1,b_2,b_3,\dots$ be the sequence of (first) differences of the $a_i$.  So $b_i=a_{i+1}-a_i$.
In general, if the sequence $b_0,b_1,b_2,\dots$ is given by $b_i=P(i)$, where $P$ is a polynomial of degree $k$, then the original sequence is given by a polynomial of degree $k+1$. This is a result in a large and useful subject called the calculus of finite differences.
It is analogous to the result that if the derivative of a function is a polynomial of degree $k\ge 0$, then the function is a polynomial of degree $k+1$. 
Thus, automatically, the sequence of the post is given by a polynomial of degree $3$. 
A: Not quite an answer but a partial explanation: the equation can be rewritten as $T(n)-T(n-1) = n(n-1)$ (note that the Fibonacci equation can't readily be 'cleared' in this fashion); we know if $p(x)$ is polynomial of degree $n$ then the finite difference $\Delta p(x) \equiv p(x+1)-p(x)$ is polynomial of degree $n-1$ (this isn't too hard to prove using the binomial theorem and makes a nice expression — you should even be able to get an explicit form for the coefficients of $\Delta p(x)$ if you try), so it makes sense to try and 'invert' this and work with the assumption that $T()$ is a third-degree polynomial for the given problem.
In fact, the so-called Bernoulli polynomials give one straightforward approach to solving this problem: they serve to 'invert' the finite-difference operator for monomials $x^n$ (in other words, if $B_{n+1}(x)$ is the order $n+1$ Bernoulli polynomial, then we have $\Delta B_{n+1}(x) = x^n$, up to a small multiplicative factor) in the same way that regular monomials do the same thing for derivatives.  Since the finite-difference operator is linear (that is, $\Delta(ap_1(x)+bp_2(x)) = a\Delta p_1(x)+b\Delta p_2(x)$), once we understand how to compute the inverse of the finite-difference operation for each individual term of a polynomial we can compute the inverse on the polynomial as a whole.
