Another way to prove it using basic linear algebra would be as follows.
First of all, the set [itex]E[/itex] of sequences [itex](u_n)[/itex] such as [itex]u_{n+2} = u_{n+1} + u_n[/itex] for all [itex]n \in \mathbf N[/itex] is a vector space over [itex]\mathbf R[/itex] (this is very easy to show).
Also, [itex]E[/itex] is of dimension 2 (intuitively, any sequence of [itex]E[/itex] is entirely defined by its two first terms, i.e. [itex]E[/itex] is isomorphic to [itex]\mathbf R^2[/itex]).
Next, by calling [itex](F_n)[/itex] our Fibonacci sequence, [itex](F_n) \in E[/itex], so [itex](F_n)[/itex] can be expressed as a linear combination of two basis vectors (here, sequences) of [itex]E[/itex].
This is the tricky part. You have to find two sequences, not proportional to each other, verifying the recurrence relation [itex]u_{n+2} = u_{n+1} + u_n[/itex], that you can clearly explicit. You can easily show that [itex](\varphi^n)[/itex] and [itex](\bar{\varphi}^n)[/itex] verify that relation, where [itex]\varphi[/itex] and [itex]\bar{\varphi}[/itex] are the solutions to the equation [itex]x^2=x+1[/itex], namely [itex]\dfrac{1+\sqrt 5}{2}[/itex] and [itex]\dfrac{1-\sqrt 5}{2}[/itex] .
Indeed, [itex]\varphi^2 = \varphi + 1[/itex] by definition, and multiplying by [itex]\varphi^n[/itex] gives [itex]\varphi^{n+2} = \varphi^{n+1} + \varphi^n[/itex] (same goes for [itex]\bar{\varphi}[/itex]).
So the sequences [itex](\varphi^n)[/itex] and [itex](\bar{\varphi}^n)[/itex] verify the recurrence relation, i.e. are vectors of [itex]E[/itex]. They aren't proportional (you can see that by checking the two first terms of each sequence), so they form a basis of [itex]E[/itex]. Thus, you can express [itex](F_n)[/itex] as a linear combination of [itex](\varphi^n)[/itex] and [itex](\bar{\varphi}^n)[/itex], in other words, there exists [itex]\lambda, \mu[/itex] reals such that for all [itex]n[/itex], [itex]F_n = \lambda \varphi^n + \mu\bar{\varphi}^n[/itex].
Evaluating at [itex]n = 0[/itex] and [itex]n = 1[/itex] gives [itex]\lambda + \mu = 0[/itex] and [itex]\lambda\varphi + \mu\bar{\varphi} = 1[/itex], which solves into [itex]\lambda = -\mu = \dfrac{1}{\varphi-\bar{\varphi}}[/itex].
Finally, [itex]F_n = \dfrac{\varphi^n - \bar{\varphi}^n}{\varphi-\bar{\varphi}}[/itex].