Proof for the formula of the $n^\text{th}$ term of a linear and homogeneous second-order recurrence

Expanding on my comments under OP's question, suppose we do not know about the characteristic polynomial, and do not have the foresight to try to find solutions in the form of a geometric progression. Then the $2^{nd}$ order recurrence can still be solved by reducing it to the known case of a $1^{st}$ order recurrence $g_n = \mu g_{n-1}$ which is easily solved by telescoping as $g_n = \mu g_{n-1}=\mu^2g_{n-2}=\dots=\mu^ng_0\,$.

Let the $2^{nd}$ order homogeneous recurrence be $f_n=a f_{n-1}-bf_{n-2}$. We can assume WLOG that $b\ne0$, otherwise it would not be a $2^{nd}$ order recurrence.

The idea is to try to rewrite the relation $f_n=a f_{n-1}-bf_{n-2}$ in the form $g_n = \mu g_{n-1}$ where $g_n=f_n-\lambda f_{n-1}$ is a linear combination between consecutive terms of the original sequence:

$$ g_n = \mu g_{n-1} \;\iff\; f_n - \lambda f_{n-1} = \mu(f_{n-1}-\lambda f_{n-2}) \;\iff\; f_n = \color{red}{(\lambda+\mu)}\,f_{n-1}- \color{green}{\lambda\mu}\, f_{n-2} $$

The last equality matches $f_n=\color{red}{a} f_{n-1}-\color{green}{b}f_{n-2}$ if we choose $\lambda,\mu$ such that $\lambda+\mu=a$ and $\lambda\mu = b$ $\;\left(\dagger\right)\;$ with $\lambda,\mu \ne 0$ because $b \ne 0$.

Then $g_n=\mu^{n-1} g_1\,$, and $f_n$ can be determined by adding the following and telescoping again:

$$ \begin{cases} f_n - \lambda f_{n-1} &= \mu^{n-1}\,g_1 && \mid \,\cdot\,1 \\ f_{n-1} - \lambda f_{n-2} &= \mu^{n-2}\,g_1 && \mid \,\cdot\, \lambda \\ f_{n-2} - \lambda f_{n-3} &= \mu^{n-3}\,g_1 && \mid \,\cdot\, \lambda^2 \\ \dots \\ f_2 - \lambda f_1 &= \mu \,g_1 && \mid \,\cdot\, \lambda^{n-2} \\ f_1 - \lambda f_0 &= \;\;\,g_1 && \mid \,\cdot\, \lambda^{n-1} \end{cases} $$

$$ \implies\;\;\;\;f_n = \lambda^n f_0 + \underbrace{\left(\lambda^{n-1}+\lambda^{n-2}\mu+\dots+\mu^{n-1}\right)}_{=\;\begin{cases}\begin{align} (\lambda^n-\mu^n)/(\lambda-\mu) && \text{if}\;\; \lambda \ne \mu \\ n\,\lambda^{n-1} && \text{if}\;\; \lambda = \mu\end{align}\end{cases}}\,\left(f_1-\lambda f_0\right) $$

$\;\left(\dagger\right)$ Note that, by Vieta's relations, these $\lambda, \mu$ are the roots of the quadratic $t^2 - at + b = 0$, which "happens" to be the characteristic polynomial associated with the original recurrence.

If you want a proof without using induction here is one. Let $r\in\mathbb{R}$ and $u_n:=r^n$, then $(u_n)$ satisfies the recursive relation $u_{n+2}=au_{n+1}+bu_n$ if and only if $r^2=ar+b$. Suppose there are 2 distinct solutions for this equation which I will denote $r_1$ and $r_2$. Then $(r_1^n)$ and $(r_2^n)$ are in E, the set of the sequences satisfying the recursive formula. They are linearly independent because $r_1\neq r_2$ thus $\dim{\rm span}((r_1^n),(r_2^n))=2$. But $\dim E\leqslant 2$ because the map $(u_n)\in E\mapsto (u_0,u_1)\in\mathbb{R}^2$ is one-to-one (you can show it by induction). We can therefore conclude that $E={\rm span}((r_1^n),(r_2^n))$ by an argument of dimension.

EDIT (other solution without linear algebra) :

Let $r_1$ be a root of $X^2-aX-b$ and $v_n:=u_nr_1^{-n}$, then you can check that $$ v_{n+2}-v_{n+1}=-\frac{b}{r_1^2}(v_{n+1}-v_n) $$ Therefore $(v_{n+1}-v_n)$ is a geometric sequence : $$ v_{n+1}-v_n=C(r_2/r_1)^n $$ where $C$ is a constant and $r_2:=-\frac{b}{r_1}$ (notice that $r_1$ and $r_2$ are the roots of $X^2-aX-b$). Summing this gives that $$ v_n=\lambda+\mu(r_2/r_1)^n $$ for some constants $\lambda,\mu$ so that $$ u_n=\lambda r_1^n+\mu r_2^n $$

Here is something I worked out a few years ago. It might be too general.

This is a totally non-original working out of the general homogeneous linear recurrence.

(Added later: the special cases of one and two term recurrences are shown explicitly at the end.)

If $a_n =\sum_{k=1}^m c_k a_{n-k} $ for $n \ge m$, find the generating function $\sum_{n=0}^{\infty} a_n x^n $ in terms of $a_{0..m-1}$ and $c_{1..m}$.

(This is undoubtedly a duplicate, but I just wanted to work this out from scratch. The algebra and rearrangements of summations is complicated enough so that, as usual, I am hoping that someone might have a simpler derivation.)

Given the linear recurrence

$a_n =\sum_{k=1}^m c_k a_{n-k} $ for $n \ge m$.

Find the generating function $A(x) =\sum_{n=0}^{\infty} a_n x^n $ in terms of $a_{0..m-1}$ and $c_{1..m}$.

My answer is

$A(x) =\dfrac{\sum_{r=0}^{m-1} x^r(a_r-\sum_{k=1}^{r} c_{k} a_{r-k})}{1-\sum_{k=1}^m c_k x^k} $.

My derivation.

Let $A_j(x) =\sum_{n=0}^{j} a_n x^n $ for $j = 0$ to $m-1$.

$\begin{array}\\ A(x) &=\sum_{n=0}^{\infty} a_n x^n\\ &=\sum_{n=0}^{m-1} a_n x^n+\sum_{n=m}^{\infty} a_n x^n\\ &=A_{m-1}(x)+\sum_{n=m}^{\infty} x^n\sum_{k=1}^m c_k a_{n-k}\\ &=A_{m-1}(x)+\sum_{k=1}^m c_k\sum_{n=m}^{\infty} x^n a_{n-k}\\ &=A_{m-1}(x)+\sum_{k=1}^m c_k x^k\sum_{n=m}^{\infty} x^{n-k} a_{n-k}\\ &=A_{m-1}(x)+\sum_{k=1}^m c_k x^k\sum_{n=m-k}^{\infty} x^{n} a_{n}\\ &=A_{m-1}(x)+\sum_{k=1}^m c_k x^k(\sum_{n=0}^{\infty} x^{n} a_{n}-\sum_{n=0}^{m-k-1} x^{n} a_{n})\\ &=A_{m-1}(x)+\sum_{k=1}^m c_k x^k(A(x)-A_{m-k-1}(x))\\ &=A_{m-1}(x)+\sum_{k=1}^m c_k x^kA(x)-\sum_{k=1}^m c_k x^kA_{m-k-1}(x)\\ \end{array} $

so $\begin{array}\\ A(x)(1-\sum_{k=1}^m c_k x^k) &=A_{m-1}(x)-\sum_{k=1}^m c_k x^kA_{m-k-1}(x)\\ &=A_{m-1}(x)-\sum_{k=0}^{m-1} c_{m-k} x^{m-k}A_{k-1}(x)\\ &=A_{m-1}(x)-\sum_{k=0}^{m-1} c_{m-k} x^{m-k}\sum_{n=0}^{k-1} a_n x^n\\ &=A_{m-1}(x)-\sum_{n=0}^{m-2}\sum_{k=n+1}^{m-1} c_{m-k} x^{m-k} a_n x^n\\ &=A_{m-1}(x)-\sum_{n=0}^{m-2}\sum_{k=n+1}^{m-1} c_{m-k} a_n x^{n+m-k}\\ &=A_{m-1}(x)-\sum_{n=0}^{m-2}\sum_{r=n+m-m+1}^{n+m-n-1} c_{m-n-m+r} a_n x^{r} \qquad r = n+m-k, k = n+m-r\\ &=A_{m-1}(x)-\sum_{n=0}^{m-2}\sum_{r=n+1}^{m-1} c_{r-n} a_n x^{r}\\ &=A_{m-1}(x)-\sum_{r=1}^{m-1}\sum_{n=0}^{r-1} c_{r-n} a_n x^{r}\\ &=\sum_{r=0}^{m-1} a_r x^r-\sum_{r=1}^{m-1}x^r\sum_{n=0}^{r-1} c_{r-n} a_n\\ &=a_0+\sum_{r=1}^{m-1} x^r(a_r-\sum_{n=0}^{r-1} c_{r-n} a_n)\\ &=\sum_{r=0}^{m-1} x^r(a_r-\sum_{n=0}^{r-1} c_{r-n} a_n)\\ &=\sum_{r=0}^{m-1} x^r(a_r-\sum_{n=1}^{r} c_{n} a_{r-n})\\ &=\sum_{r=0}^{m-1} x^r(a_r-\sum_{k=1}^{r} c_{k} a_{r-k})\\ \end{array} $

Therefore

$A(x) =\dfrac{\sum_{r=0}^{m-1} x^r(a_r-\sum_{k=1}^{r} c_{k} a_{r-k})}{1-\sum_{k=1}^m c_k x^k} $.

Once you have this, you can use partial fractions to decompose the denominator in terms of the roots.

(added later)

For $m=1$, when $a_n =c_1a_{n-1} $ for $n \ge 1 $, the initial term is $a_0$, and the expressions in braces ("{...}") give the values of the indices, $$A(x) =\dfrac{a_0 }{1-c_1x} =\dfrac{\{r=0\}a_0 }{1-c_1x} $$

For $m=2$, when $a_n =c_1a_{n-1}+c_2a_{n-2} $ for $n \ge 2 $, the initial terms are $a_0$ and $a_1$, and the expressions in braces ("{...}") give the values of the indices, $$A(x) =\dfrac{a_0 +x(a_1-c_1a_{0})}{1-c_1x-c_2x^2} =\dfrac{\{r=0\}a_0 +\{r=1\}x(a_1-\{k=1\}c_1a_{0})}{1-c_1x-c_2x^2} $$