You have already answered your own question, essentially. The space $E$ of sequences $\{a_n\}_{n\geq 0}$ fulfilling $a_{n+2}=a_{n+1}+a_n$ is a vector space with dimension $2$; any sequence belonging to $E$ is fixed by its initial values $a_0$ and $a_1$. The initial values of the Fibonacci sequence are $0,1$; the initial values of the Lucas sequence are $2,1$; since $(0,1)$ and $(2,1)$ are linearly independent each sequence in $E$ can be represented through $a_n = f\cdot F_n + l\cdot L_n$ for some constants $f,l$. For instance, this applies to the shifted Fibonacci and Lucas sequences: $$ L_{n+1} = \frac{5}{2}F_n+\frac{1}{2}L_n,\qquad F_{n+1}=\frac{1}{2}F_n+\frac{1}{2}L_n $$ and you may also take $\{F_n\}_{n\geq 0},\{F_{n+1}\}_{n\geq 0}$ or $\{L_n\}_{n\geq 0},\{L_{n+1}\}_{n\geq 0}$ as a base for $E$.
The canonical choice should have been $(1,0),(0,1)$ associated to $\{F_{n-1}\}_{n\geq 0},\{F_{n}\}_{n\geq 0}$, but there are arithmetic reasons for picking $\{F_n\}_{n\geq 0}$ and $\{L_n\}_{n\geq 0}$ as the "mostly peculiar" elements of $E$ (see the comments below).
A sequence with a moderate growth can be associated to an analytic function via $$ \{a_n\}_{n\geq 0}\quad \mapsto\quad g(x)=\sum_{n\geq 0}a_n x^n $$ (the RHS is known as the ordinary generating function (OGF) of the sequence) and this gives a tight relation between linear recurrent sequences, linear differential equations and the elements of $M^n$ with $M$ being a fixed matrix. You guessed it correctly: they all appear to be the same problem because they actually are the same problem.
Here it is a brief outline of the power of the previous association $\text{sequence}\mapsto\text{OGF}$, relying on the fact that we may go in the opposite direction by Cauchy's integral formula or something equivalent to it. Since in the Fibonacci sequence each term with $n\geq 2$ is given by the sum of the previous two terms, when we multiply the $\text{OGF}$ by $1-x-x^2$ there has to occur a massive cancellation, leaving out at most two monomials. Indeed, by interpolation: $$ \sum_{n\geq 0}F_n x^n = \frac{x}{1-x-x^2}.$$ This is beautiful, since we may factor the polynomial $1-x-x^2$, apply a partial fraction decomposition to the RHS and recall that $$ \frac{1}{1-\alpha x} = \sum_{n\geq 0} \alpha^n x^n $$ holds for any $\alpha$ such that $|\alpha|<1$. By the unicity of the Maclaurin series of the $\text{OGF}$, the explicit formula $$ F_n = \frac{1}{\sqrt{5}}\left[\left(\frac{1+\sqrt{5}}{2}\right)^n - \left(\frac{1-\sqrt{5}}{2}\right)^n\right] $$ is proved. We may apply the same technique to the sequence of Lucas numbers: the only different step is the interpolation step, leading to $$ L_n = \left[\left(\frac{1+\sqrt{5}}{2}\right)^n + \left(\frac{1-\sqrt{5}}{2}\right)^n\right]=\text{Tr}\begin{pmatrix} 0 & 1 \\ 1 & 1 \end{pmatrix}^n. $$ If we replace the constants $\varphi,1-\varphi$ with the constants $e^i, e^{-i}$ we have that $F_n$ is "essentially a (hyperbolic) sine" and $L_n$ is "essentially a (hyperbolic) cosine". Rigorously, $$\boxed{ F_n = \frac{2}{\sqrt{5}}\,\sinh(n\log\varphi),\qquad L_n = 2\cosh(n\log\varphi)\,}$$ and this also allows to define $F_n$ or $L_n$ for non-integer values of $n$.