How to *really* solve a non-homogeneous recurrence First let me state that I am not asking about the usual procedure for finding a trial solution to a non-homogeneous recurrence.  I have been doing this for many years and can solve all the basic types, but I am looking for some deeper insight.
Here are a few examples to serve as a basis for the discussion.


*

*For $a_n-a_{n-1}-6a_{n-2}=n^2$ we guess $a_n=cn^2+dn+e$.  Why not just $cn^2$?  If you want to say "because $(n-1)^2$ contains lower order terms" please read the rest of the question before posting that answer.

*For $a_n-a_{n-1}-6a_{n-2}=2^n$ we guess $a_n=c2^n$: no further terms as in the previous example.

*For $a_n-a_{n-1}-6a_{n-2}=3^n$ we guess $a_n=cn3^n$.  Why not just $c3^n$?  (And again, I do not want the answer "because $3^n$ satisfies the homogeneous recurrence", I know that already.)  Why is it not necessary to include a term $dn^23^n$?

*For $a_n-a_{n-1}-6a_{n-2}=2^n+3^n$ we guess $a_n=c2^n+dn3^n$.  How do we really know that $c2^n$ is OK but $d3^n$ is not?

*For $a_n-6a_{n-1}+9a_{n-2}=n3^n$ we guess $a_n=cn^33^n+dn^23^n$.  As above, how do we know in advance that we will not need a term $en^43^n$?
My thoughts on this are very vague, any insight (perhaps even proofs) would be greatly appreciated.
I feel that the answer must have something to do with linear independence in the vector space $V$ of all (let's say real) sequences $\{a_0,a_1,\ldots\}$.
In case 5, for example, if we start with the homogeneous solutions, we have to continue the set
$$\{3^n,\,n3^n,\ldots\}$$
until we get four independent sequences - but why four?
In case 4, I assume that we treat the two summands separately because the sequences $\{2^n\}$ and $\{3^n\}$ are linearly independent.
On the other hand, $\{n\}$ and $\{n^2\}$ are also linearly independent and this would be different.  For this reason, and also to account for example 1, I suspect that we also need to consider the finite difference operator
$$S:\{a_0,a_1,\ldots\}\mapsto\{a_1-a_0,a_2-a_1,\ldots\}\ ,$$
which is a linear transformation on $V$, and perhaps to consider whether $a_n,\,S(a_n),\,S^2(a_n),\ldots$ are linearly independent.
What I would like to see is a theorem of the following shape.

If the $k$th order linear recurrence with constant coefficients
  $$a_n+\cdots+c_ka_{n-k}=0$$
  has solutions ${\rm span}({\bf h}_1,\ldots,{\bf h}_k)$, then the recurrence
  $$a_n+\cdots+c_ka_{n-k}=f(n)$$
  has solutions of the form $a_n={}????$

Any ideas will be read with interest! - but once again, please do not tell me how to solve recurrences!  I know how, I am trying to understand why.
 A: My suggestion is to think of the undetermined coefficient method for linear non-homogeneous ODE, which can be found in any elementary ODE book. Let me illustrate. Let say we want to find a particular solution of
$$
ay''+by'+cy=f(x),
$$
where $a,b,c$ are constants. If $f(x)=x^3$, you will need to plug in a polynomial of degree 3 (try it and see why you need the lower degree terms). Now, if $f(x)=e^{ax}$, then you will need to plug in $y=Ke^{ax}$. I believe it is obvious that you do not need any other terms. Notice that if $f(x)=2^x$, it is a particular case of an exponential since $f(x)=2^x=e^{\ln{(2)} x}$. Now, the only exceptions to those rules if if $f(x)$ is a solution to the homogeneous problem. For example in the case of
$$
y''-y=e^x.
$$
In this case, $y=ke^x$ does not work. You need to pre-multiply by $x$ to try $y=kxe^x$. In general, if $f(x)$ is a solution of the homogeneous problem, you need to pre-multiply the ansatz's I suggested above by $x$. All these results are formulated in theorems in any decent elementary ODE book (for example, the one written by Nagle and Saff). But I think these rules can be understood quite easily by working out simple examples. Note that in relation to your example $n3^n$, I can add the rule that if $f(x)=x^2 e^{ax}$ and that $e^{ax}$ is not a solution of the homogeneous problem, you need to try a polynomial of degree $3$ times $e^{ax}$. Note that if you try to add a term like $x^4e^{ax}$ to your solution, then this term will automatically produce a term proportional to $x^4e^{ax}$, which is unwanted.
If $e^{ax}$ was a solution, then you would need to pre-multiply by $x$. You would try $y=x(bx^3+cx^2+dx+f)e^{ax}$. Again, all this are in the theorems. Now, if $a$ was a double root of the characteristic equation, you would need to pre-multiply by $x^2$ instead.
For example
$$
y''-6y'+9y=x^2e^{3x}.
$$
You would try $y=x^2(ax^2+bx+c)e^{3x}.
$
Note that in the case of $n3^n$ you mention, the characteristic polynomial has a double root $3$, thus the homogeneous problem has $3^n$ and $n3^n$ as solutions, that is why you need  to pre-multiply by $n^2$ and try $n^2(an+b)3^n$,
All what I said I believe extends to your difference equations. Your difference equations are ``like'' second order differential equation with inhomogeneities made of exponentials and polynomials. So, the rules for the undetermined coefficient method do extend.  
A: Here's the linear algebra connection (loosely). Look at these three matrices: 
$$
M_1 = \pmatrix{2 & 0 \\ 0 & 3}\\
M_2 = \pmatrix{3 & 0 \\ 0 & 3}\\
M_3 = \pmatrix{3 & 1 \\ 0 & 3}\\
$$
The first has two distinct eigenvalues, so you know there's an eigenvector for each one. Nice and simple. 
$M_2$ has the same eigenvalue twice...but you can still pick two distinct eigenvectors (say $e_1$ and $e_2$). 
$M_3$ has the same eigenvalue twice, but only one eigenvector (namely $e_1$) corresponding to that eigenvector. There is, however, a "generalized eigenvector" $w$ with the property that $M_3 w$ is an eigenvector, i.e., $w$ itself isn't an eigenvector, but once you apply the transformation once, $w$ "falls into" the real eigenspace. 
If you build up a larger matrix, you can get things like
$$
\pmatrix{
2 & 1 & & & & \\
  & 2 & 1 & & & \\
   & & 2 & 1 & & \\
   & & & 2 &  & & \\
   & & & & 5 & 1 &\\
   & & & &  & 5 & &\\
   & & & &  &  & 3&\\
}
$$
where $2$ has one true eigenvector, and 3 more "generalized" eigenvectors: one transforms immediately to an eigenvector, the next, transformed twice, becomes an evec; the last, transformed 3 times, becomes an evec. 
At the same time, $5$ has one true evec, and one generalized evec, and $3$ has only a true evec. 
What's this have to do with differential equations? (Like @Stephane, I'm going to do differential rather than difference equations, because it's easier to write.)  
Well, let's first think about equations whose solutions are (possibly infinite) polynomials. There's a nice basis for polynomials (1, x, x^2, ...). What does the "matrix" for the derivative operator look like in this basis? Like this:
$$
D = \pmatrix{0 & 0 & 0 & \ldots \\
1 & 0 & 0 & \ldots \\
  & 2 & 0 & \ldots \\
  &   & 3 & \ldots \\
\vdots}
$$
Now to use the analysis that I presented before, you'd need to rewrite this in Jordan normal form, but you can see right away that there's an eigenvector for $0$ (namely $p(x) = 1$), and a generalized eigenvector for $0$ (namely $p(x) = x$) and ...lots more of those. 
Look at some ODE like 
p'' - 2 p' + p = 2x + 1
The matrix for the left hand side will be 
$$
D^2 -2D + I
$$
which may have various eigenvectors for $0$ ("solutions of the homogeneous equation") and may also have generalized eigenvectors for $0$. 
... and my willingness/ability to elaborate here has petered out. The key thing is that when you're solving ODEs, you're looking at eigenvectors for linear maps (possibly on large spaces of functions!), and something very like the Jordan Normal Form is the relevant representation of those operators to understand which ones might arise in your "general" solution, etc. 
A: $$
a_{k} - a_{k-1} - 6 a_{k-2} = k^2  \tag{1}\\
$$
$$
a_{k+1}- a_{k} - 6 a_{k-1} = (k+1)^2 \tag{2}\\
$$
denote $\Delta a_{k+1} = a_{k+1} -a_k$ and subtract (1) from (2) to get 
$$
\Delta a_{k+1} - \Delta a_{k} - 6 \Delta a_{k-1} = 2k -1 \tag{3}\\
$$
Now 
$$
\Delta a_{k+2} - \Delta a_{k+1} - 6 \Delta a_{k} = 2(k+1) -1 \tag{4}\\
$$
denote $b_{k+2} = \Delta a_{k+2} - \Delta a_{k+1}$ and subtract (3) from (4) to get
$$
\Delta b_{k+2} - \Delta b_{k+1}-6 \Delta b_k = 2
$$
Thus you've got rid of the nonhomogenous term on the RHS. Now rewrite the expression as ($V_{k+2} = \Delta b_{k+2} - \Delta b_{k+1}) $
$$
V_{k+2} = 2  + 6 \Delta b_{k}
$$
and use telescoping sums on both sides (sum on $k$). I know this is slower than the approach you've suggested, but at least it doesn't involve guessing.
A: Ask, and ye shall receive.
Shift indices to write the recurrence as:
$$
a_{n + 2} - a_{n + 1} - 6 a_n = n^2 + 4 n + 4
$$
Define the generating function $A(z) = \sum_{n \ge 0} a_n z^n$, multiply your recurrence by $z^n$ and sum over $n \ge 0$. Recognize some sums to get:
$$
\frac{A(z) - a_0 - a_1 z}{z^2}
  - \frac{A(z) - a_0}{z}
  - 6 A(z)
  = \sum_{n \ge 0} n^2 z^n + 4 \sum_{n \ge 0} n z^n + 4 \sum_{n \ge 0} z^n
$$
The tricky part here are the sums on the left hand side. But you see that multiplying the coefficient of a series $\sum_{n \ge 0} c_n z^n$ by $n$ results from differentiating and multiplying by $z$ to restore the power "lost" to the derivative. So:
$\begin{align}
\sum_{n \ge 0} z^n
  &= \frac{1}{1 - z} \\
\sum_{n \ge 0} n z^n
  &= z \frac{\mathrm{d}}{\mathrm{d} z} \frac{1}{1 - z} \\
  &= \frac{z}{(1 - z)^2} \\
\sum_{n \ge 0} n^2 z^n
  &= z \frac{\mathrm{d}}{\mathrm{d} z} \frac{z}{(1 - z)^2} \\
  &= \frac{z}{(1 - z)^2} \\
  &= \frac{z + z^2}{(1 - z)^3}                
\end{align}$
Plug all this into the above, and you'll get an algebraic equation for $A(z)$, which gives $A(z)$ as a rational function in $z$ in this case. Split into partial fractions, and remember that:
$$
(1 - r z)^{-k}
  = \sum_{n \ge 0} (-1)^n \binom{-k}{n} r^n z^n
  = \sum_{n \ge 0} \binom{n + k - 1}{k - 1} r^n z^n
$$
This allows you to "read off" the coefficients.
The real fun starts when you have complex roots, which appear in conjugate pairs, call one set of them $\alpha$ and $\overline{\alpha}$, with corresponding terms of the partial fraction expansion (again $\beta$ and $\overline{\beta}$ are complex conjugates):
$$
\frac{\beta}{1 - \alpha z}
  + \frac{\overline{\beta}}{1 - \overline{\alpha} z}
$$
This gives rise to a pair of terms in the solution:
$$
\beta \cdot \alpha^n + \overline{\beta} \cdot \overline{\alpha}^n
$$
This are seen to be the sum of complex conjugates, i.e., this is twice the real part of any of them:
$$
2 \Re \left( \beta \cdot \alpha^n \right)
$$
Write the numbers is polar form, i.e., $\alpha = \lvert \alpha \rvert \cdot \exp(\arg(\alpha) \mathrm{i})$, which gives:
$$
2 \Re \left(
        \lvert \beta \rvert \exp( \arg(\beta) \mathrm{i} )
          \cdot \lvert \alpha \rvert^n \exp( \arg(\alpha) n \mathrm{i})
        \right)
  = 2 \lvert \beta \lvert \cdot \lvert \alpha \rvert^n
        \cdot \cos(\arg(\alpha) n + \arg(\beta))
$$
Repeat complex roots are handled in a similar way, giving terms with binomial coefficients as above.
Much more details are given in Wilf's "generatingfunctionology", at the link you'll find the next to last edition in PDF for free.
