Assuming $b_0\ne0$, a way to proceed is to introduce $A(x)=\sum\limits_{n=0}^{+\infty}a_nx^n$ and $B(x)=\sum\limits_{n=0}^kb_nx^n$ and to note that $A(x)B(x)$ has no power of $x$ of degree $\geqslant k$. Thus $A(x)=C(x)/B(x)$, where $C(x)$ is a polynomial of degree at most $k-1$ depending on $(b_n)_{0\leqslant n\leqslant k-1}$ and $(a_n)_{0\leqslant n\leqslant k-1}$.
The rest is as you describe it: decompose the rational function $C(x)/B(x)$ and collect the coefficients $a_n$ for $n\geqslant k$. For every root $1/r$ of $B$ with multiplicity $m$, the decomposition of $C(x)/B(x)$ includes terms $1/(1-rx)^s$ for $s\leqslant m$, and one uses the fact that
$$
\frac1{(1-rx)^s}=\sum_{n=0}^{+\infty}{n+s\choose s}r^nx^n,
$$
and that ${n+s\choose s}$ is a polynomial in $n$ with degree $s$ to conclude.
A nice proof is given by using generating functions. Say your recurrence is:
$\begin{align*}
c_k a_{n + k} + c_{k - 1} a_{n + k - 1} + \dotsb + c_0 a_n = 0
\end{align*}$
with given values for $a_0, \dotsc, a_{k - 1}$. Define the generating function $A(z) = \sum_{n \ge 1} a_n z^n$, multiply your recurrence by $z^n$ and add up over $n \ge 0$. Recognizing some sums you have:
$\begin{align*}
c_k \frac{A(z) - a_0 - a_1 z - \dotsb - a_{k - 1} z^{k - 1}}{z^k}
+ c_{k - 1} \frac{A(z) - a_0 - a_1 z - \dotsb - a_{k - 2} z^{k - 2}}{z^{k - 1}}
+ \dotsb
+ c_0 A(z)
&= 0
\end{align*}$
Multiply through by $z^k$ and solve for $A(z)$ to get:
$\begin{align*}
A(z)
&= \frac{p(z)}{c_k z^k + c_{k - 1} z^{k - 1} + \dotsb + c_0}
\end{align*}$
for some polynomial $p(z)$ of degree at most $k - 1$ that depends on the initial values. Note that the denominator is the characteristic equation of the recurrence.
From calculus you know that this can be written as partial fractions:
$\begin{align*}
A(z)
&= \sum_r \frac{A_r}{(1 - \alpha_r z)^{m_r}}
\end{align*}$
where $1 / \alpha_r$ is a zero of the characteristic equation of multiplicity $m_r$. Now, by the extended binomial theorem:
$\begin{align*}
(1 - \alpha z)^{-m}
&= \sum_{n \ge 0} (-1)^n \binom{-m}{n} (\alpha z)^n \\
&= \sum_{n \ge 0} \binom{m + n - 1}{m - 1} \alpha^n z^n
\end{align*}$
But the binomial coefficient $\binom{m + n - 1}{m - 1}$ is just a polynomial in $n$ of degree $m - 1$. Thus the solution is a sum of terms of the form $\alpha^n$ multiplied by a polynomial in $n$ of degree at most $m - 1$ if $\alpha$ is a zero of multiplicity $m$.
If any zero turns out complex, there will be a pair of complex conjugate zeroes. Writing out the complex numbers in polar form $\alpha = \rho \exp(\omega i)$ you get $\alpha^n = \rho^n \exp(n \omega i)$, giving rise to terms of the form $\rho^n \cos(\omega n)$ (they will turn out in conjugate pairs, the imaginary parts cancel out).
Best Answer
The usual convention is $f_0=0, f_1=f_2=1$. Using your equations $$c_1+c_2=0\\c_1r_1+c_2r_2=1\\ c_1\frac{1+\sqrt 5}2+c_2\frac{1-\sqrt 5}2=1\\ \frac {\sqrt 5}2(c_1-c_2)=1\\ c_1=\frac 1{\sqrt 5}\\c_2=-\frac 1{\sqrt 5}$$