Why this method of solving recurrence relation works? Could anyone explain why we can solve recurrence relations by finding the soltuion of its characteristic equation? I'm talking about the method presented here. Is the proof of the method validity so complex that it has been omitted in all texts I've come across?
Does it always work? Example:
$a_n = c_{1}a_{n-1} + c_{2}a_{n-2}+...+c_{k}a_{n-k}$
Then $r^{n}$ is a solution to above equation if and only if $r$ is a solution of $r^{k}-c_{1}r^{k-1} - c_{2}^r{k-2}-...-c_k = 0$.
 A: No, it’s pretty straightforward. The case $n=2$ is proved in most of the elementary textbooks that I’ve seen.
Suppose that $a_n=r^n$ is a solution to the recurrence
$$a_n=\sum_{i=1}^kc_ia_{n-i}\;.\tag{1}$$
Then 
$$r^n=\sum_{i=1}^kc_ir^{n-i}=r^{n-k}\sum_{i=1}^kc_ir^{k-i}\tag{2}$$
for $n\ge k$. In particular, setting $n=k$ we have
$$r^k=\sum_{i=1}^kc_ir^{k-i}\;,\tag{3}$$
i.e.,
$$r^k-c_1r^{k-1}-c_2r^{k-2}-\ldots-c_k=0\;.\tag{4}$$
Conversely, if $r$ satisfies $(4)$, then it satisfies the rewritten version $(3)$. For any $n\ge k$ we can multiply $(3)$ by $r^{n-k}$ to find that $r$ satisfies $(2)$ for every $n\ge k$, and it follows immediately that $a_n=r^n$ is a solution to $(1)$.
A: Indeed $r^n$ is a solution to this recursion iff it is a root the characteristic polynomial.  As Brian M. Scott did in his answer you can just plug it into the recurrence and solve.  Moreover if $r$ is a $k$-tuple root of the characteristic polynomial then $r^n, nr^n, n^2r^n,... n^{k-1}r^n$ are all solutions, this can be seen by repeatedly differentiating $x^n$ times the characteristic polynomial and remembering that a $k$-tuple root is also the root of the first $k-1$ derivatives as well.
What's more is that any solution to such a recurrence is a linear combination of the sequences found this way. The proofs I know of this all use some amount of linear algebra, so it's not surprising to me that they are often omitted in elementary explanations like the one you linked to.  Depending on how comfortable you are with linear algebra I can elaborate on this more if you'd like.
A: An idea why it works is given by generating functions. Take the recurrence:
$\begin{equation*}
c_r a_{n + r} + c_{r - 1} a_{n + r - 1} + \dotsb + c_c a_n
  = f(n)
\end{equation*}$
Define:
$\begin{equation*}
  A(z)
    = \sum_{n \ge 0} a_n z^n
\end{equation*}$
Multiply the recurrence by $z^n$, sum over $n \ge 0$ and recognize the resulting sums:
$\begin{align*}
  c_r \sum_{n \ge 0} a_{n + r} z^n
    + \dotsb
    + c_0 \sum_{n \ge 0} a_n z^n
   &= \sum_{n \ge 0} f(n) z^n \\
  c_r 
   \frac{A(z) - a_0 - a_1 z - \dotsb - a_{r - 1} z^{r - 1}}
        {z^r}
     + \dotsb
     + c_0 A(z)
    &= F(z)
\end{align*}$
Multiply through by $z^r$, you get an expression of the form:
$\begin{equation*}
   (c_r + c_{r - 1} z + \dotsb c_0 z^r) A(z)
     = z^r F(z) + p(z)
\end{equation*}$
where $p(z)$ is a polynomial depending on the initial values $a_0, \dotsc, a_{r - 1}$. If you solve for A(z), assuming $F(z)$ is a rational function, you get a rational function in $z$, of which you want the coefficients of the Maclaurin series (those are your $a_n$). The easy way to get them is to split $A(z)$ into partial fractions (don't be shy, factor down to complex factors!). You get terms of the form:
$\begin{align*}
  \frac{1}{(1 - \alpha z)^m}
    &= \sum_{k \ge 0} (-1)^k \binom{-m}{k} \alpha^k z^k \\
    &= \sum_{k \ge 0} \binom{k + m - 1}{m - 1} \alpha^k z^k 
\end{align*}$
Here $1 / \alpha$ is a zero of
$c_r + c_{r - 1} z + \dotsb + c_0 z^n$,
i.e., $\alpha$ is a zero of
$c_r \rho^r + c_{r - 1} \rho^{r - 1} + \dotsb + c_0$.
Also note that the binomial coefficient $\binom{k + m - 1}{m - 1}$ is a polynomial in $k$ of degree $m - 1$ (thus the formula for case of multiple zeros in the characteristic equation). If some factors of the characteristic equation repeat in the denominator of $F(z)$, they will add multiplicities to the respective factors, adding powers of $n$.
In case you get complex zeros, they show up as complex conjugates, say $\alpha, \bar{\alpha}$. Then you have in the partial fraction split:
$\begin{equation*}
  \frac{c}{1 - \alpha z}
    + \frac{\bar{c}}{1 - \bar{\alpha} z}
\end{equation*}$
This gives rise to terms in the expansion:
$\begin{equation*}
  c \alpha^n + \bar{c} \bar{\alpha}^n
\end{equation*}$
They are complex conjugates, so their imaginary parts cancel out. If you write in polar coordinates, by Euler's equation:
$\begin{align*}
  \alpha
    &= \rho (\cos \theta + i \sin \theta) \\
  \alpha^n
    &= \rho^n (\cos n \theta + i \sin n \theta)
\end{align*}$
If you do the whole dance of writing out and simplifying the respective term pair, you see it is real and can be written:
$\begin{equation*}
  c' \rho^n \cos (n \theta + \phi)
\end{equation*}$
where $c', \phi$ depend on $c$. Higher multiplicities give rise to to polynomials in $n$ multiplying this, just as above.
