What are some mathematically interesting computations involving matrices? I am helping designing a course module that teaches basic python programming to applied math undergraduates. As a result, I'm looking for examples of mathematically interesting computations involving matrices. 
Preferably these examples would be easy to implement in a computer program. 
For instance, suppose
$$\begin{eqnarray}
F_0&=&0\\
F_1&=&1\\
F_{n+1}&=&F_n+F_{n-1},
\end{eqnarray}$$
so that $F_n$ is the $n^{th}$ term in the Fibonacci sequence. If we set
$$A=\begin{pmatrix}
1 & 1 \\ 1 & 0 
\end{pmatrix}$$
we see that
$$A^1=\begin{pmatrix}
1 & 1 \\ 1 & 0
\end{pmatrix} =
\begin{pmatrix}
F_2 & F_1 \\ F_1 & F_0
\end{pmatrix},$$
and it can be shown that
$$
A^n =
\begin{pmatrix}
F_{n+1} & F_{n} \\ F_{n} & F_{n-1}
\end{pmatrix}.$$
This example is "interesting" in that it provides a novel way to compute the Fibonacci sequence. It is also relatively easy to implement a simple program to verify the above.
Other examples like this will be much appreciated. 

Rather than give a new example, I'd like to present a way to expand on your own example. You mention that 
$$A=\begin{pmatrix}
1 & 1 \\ 1 & 0
\end{pmatrix} =
\begin{pmatrix}
F_2 & F_1 \\ F_1 & F_0
\end{pmatrix},\quad
A^n =
\begin{pmatrix}
F_{n+1} & F_{n} \\ F_{n} & F_{n-1}
\end{pmatrix}.$$
Though, as @Henrik mentions, the matrix multiplication simply consists of doing the "usual" addition of Fibonacci numbers.
However, when you see a matrix power, e.g. $A^n$, you think of diagonalisation to efficiently compute this matrix power.
Any decent mathematical software should be able to compute for you that
$$A=\begin{pmatrix}
1 & 1 \\ 1 & 0
\end{pmatrix} = \begin{pmatrix}
-\frac{1}{\varphi} & \varphi \\ 1 & 1
\end{pmatrix}\begin{pmatrix}
-\frac{1}{\varphi} & 0 \\ 0 & \varphi
\end{pmatrix}\begin{pmatrix}
-\frac{1}{\sqrt{5}} & \frac{\varphi}{\sqrt{5}} \\ \frac{1}{\sqrt{5}} & \frac{1}{\varphi\sqrt{5}}
\end{pmatrix} = P\Lambda P^{-1},$$
where $\varphi = \frac{1+\sqrt{5}}{2}$ is the golden ratio. Therefore,
$$A^n = P\Lambda^n P^{-1} = \begin{pmatrix}
-\frac{1}{\varphi} & \varphi \\ 1 & 1
\end{pmatrix}\begin{pmatrix}
\left(-\varphi\right)^{-n} & 0 \\ 0 & \varphi^n
\end{pmatrix}\begin{pmatrix}
-\frac{1}{\sqrt{5}} & \frac{\varphi}{\sqrt{5}} \\ \frac{1}{\sqrt{5}} & \frac{1}{\varphi\sqrt{5}}
\end{pmatrix}.$$
Working out the product of these three matrices is really easy, so we get:
$$\begin{pmatrix}
F_{n+1} & F_{n} \\ F_{n} & F_{n-1}
\end{pmatrix} = A^n = \begin{pmatrix}
\frac{\varphi^{n+1} - \left(-\varphi\right)^{-(n+1)}}{\sqrt{5}} & \frac{\varphi^{n} - \left(-\varphi\right)^{-n}}{\sqrt{5}} \\ \frac{\varphi^{n} - \left(-\varphi\right)^{-n}}{\sqrt{5}} & \frac{\varphi^{n-1} - \left(-\varphi\right)^{-(n-1)}}{\sqrt{5}}
\end{pmatrix}.$$
In other words, through the use of diagonalisation we have proven that
$$F_n = \frac{\varphi^{n} - \left(-\varphi\right)^{-n}}{\sqrt{5}}.$$
