This is something I have never found a convincing answer to; maybe I haven't looked in the right places. When solving a linear difference equation we usually put $u_n=r^n$, solve the resulting polynomial, and then use these to form a complete solution. This is how to get $F_n=\frac{1}{\sqrt{5}}(\phi^n-(-\phi)^{-n})$ from $F_{n+1}=F_n+F_{n-1}$ and $F_0=0, F_1=1$ for example.
I can see that this works, but I don't understand where it comes from or why it does indeed give the full solution. I have read things about eigenfunctions and other things but they didn't explain the mechanics behind this very clearly. Could someone help me understand this?
This form of the the solution comes out rather naturally if one works in terms of generating functions. This reduces the linear difference equation to a matter of careful algebra, in rather a similar way as the Laplace transform does for linear differential equations.
The Fibonacci sequence provides a nice example of this. We have $F_{n+2}=F_{n+1}+F_n$ for $n\geq 0$ and $F_0=F_1=1$ as boundary values. We may introduce the generating function / formal power series $\mathcal{F}(x)=\sum_{n=0}^\infty F_n x^n$ where $x$ is a formal variable. Then
\begin{align}\mathcal{F}(x) &=1+x+\sum_{n=0}^\infty F_{n+2} x^{n+2}\\ &=1+x+\sum_{n=0}^\infty F_{n+1} x^{n+2}+\sum_{n=0}^\infty F_{n} x^{n+2}\\ &=1+x+x(\mathcal{F}(x)-1)+x^2 \mathcal{F}(x)\\&=1+(x+x^2)\mathcal{F}(x)\hspace{2.5 cm}\implies \mathcal{F}(x)=\frac{1}{1-x-x^2} \end{align}
One may verify that series expansion of $\mathcal{F}(x)$ produces the Fibonnaci sequence as coefficients. But we can also express this generating function via partial fractions as $$\mathcal{F}(x)=\dfrac{A_+}{1-r_+ x}+\dfrac{A_-}{1-r_- x}$$ where $A_\pm$ are appropriate coefficients and $r_\pm$ the roots of $1-x-x^2=0$ (i.e. the characteristic equation!) Expanding these as geometric series, we find coefficients $F_n = A_+ (r_+)^n+A_- (r_-)^n$. This is precisely the form we had expected.