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. 

Finite Difference Method 
IMO pure math problems are often not that attractive for undergrads and pupils. People can better relate to applied problems, even if abstracted a little, than to completely abstract math problems.
So, if you want to show them something more applied here is an example using differential equations. For the 1D case you most likely don't need anything, for the 2D case you need the concept of partial differentiation. 

The following equation e.g. models the behaviour (deformation $u$) of the membrane of a drum (2D), if you apply pressure $f$ on it. The same for 1D, which is easier to explain, since partial derivatives might not be known.
        | f
        v
 |--------------|       ⇒          |--         --|
 0              1                     --      --
                                        ------

For 2D: Let $Ω=[0,1]²$. Solve: 
$$Δu(x,y) = f(x,y) \text{ in } Ω, \qquad u(x,y)=0 \text{ on } ∂Ω.$$
For 1D:  Let $I=[0,1]$. Solve: 
$$-u''(x)=f(x) \qquad \text{ in } I, \qquad u(0)=u(1)=0.$$

Let's first look at the 1D case, since it is easier to deduce the 2D case from it afterwards. They should be able to understand the idea of 1D FDM immediately. 
The idea of "F"DM is to cover the interval $I$ with $n+2$ points, so called nodal points: 
    |-----x------x------x------x-----|
    0                                1

   x_0   x_1    x_2    x_3    x_4    x_{n+1}

          <--h-->

So we now have $$0=x_0<x_1<…<x_{n+1}=1.$$
These points should be equidistant, meaning: $$x_{i+1}-x_i = h, \qquad i=0,…,n.$$
Now since the differential equation holds for all $x∈I$ it also holds for all $x_i$:
$$-u''(x_i)=f(x_i) \qquad ∀i=1,…,n; \qquad u(0)=u(1)=0.$$
The next step of F"D"M is to approximate the derivative with the difference quotient: 
$$u''(x) = \frac{u(x-Δx)-2u(x)+u(x+Δx)}{Δx^2}$$
Hence:
$$-u''(x_i) = \frac{-u(x_{i-1})+2u(x_i)-u(x_{i+1})}{h^2}, \qquad i=1,…,n$$
And we get $n$ equations:
$$ \frac{-u(x_{i-1})+2u(x_i)-u(x_{i+1})}{h^2} = f(x_i), \qquad i=1,…,n$$
and the boundary condition is:
\begin{align*} 
 u(x_0)&=0, \\ u(x_{n+1}) &= 0
\end{align*}
Written as a matrix system it is: 
$$\frac{1}{h^2}\begin{pmatrix}
1 &0     \\
-1 &2 & -1 &   \\
 & -1 &2 & -1 &    \\
 & & \ddots &\ddots & \ddots &  \\
 & &  & -1 & 2 & -1 \\
 & &  &  & & 0  & 1 \\
\end{pmatrix}
\begin{pmatrix}u_0 \\ u_1 \\ \vdots \\ \vdots \\ u_n \\ u_{n+1} \end{pmatrix}
= 
\begin{pmatrix}f_0 \\ f_1 \\ \vdots \\ \vdots \\ f_n \\ f_{n+1} \end{pmatrix}$$
using the notation $u_i = u(x_i)$, and setting $f_0=f_{n+1}=0$.
The inner part of that matrix is symmetric and positive definite, so everything you could wish. 

Edit: It might be that you need to multiply $f$ by $-1$ to get the deformation $u$ in the correct direction. Usually one looks at $-Δu=\pm f$, because of two reasons: First the matrix given by $-Δ$ is positive definite, and second the heat equation reads $∂_tu-Δu$. 
Edit2: I noticed, that symmetric positive definite needs to be explained more. The equations of the boundary values don't have to be solved. They can be transferred to the rhs. Therefore, you actually solve the following $n×n$-system: 
$$\frac{1}{h^2}\begin{pmatrix}
2 & -1 &   \\
 -1 &2 & -1 &    \\
 & \ddots &\ddots & \ddots &  \\
 &  & -1 & 2 & -1 \\
 &  &   & -1 & 2 
\end{pmatrix}
\begin{pmatrix} u_1 \\ u_2 \\ \vdots\\ u_{n-1} \\ u_n  \end{pmatrix}
= 
\begin{pmatrix} f_1 + h^2f_0 \\ f_2 \\ \vdots \\ f_{n-1}\\ f_n + h^2f_{n+1} \end{pmatrix}$$
And this matrix is symmetric, positive definite by weak row sum criterion.

In 2D you need a grid of points, e.g. $m×m$ points.
 If you sort them lexicographically (row-wise)
        |           |           |
---  (i+m-1)------(i+m)------(i+m+1) ---
        |           |           |
        |           |           |
---   (i-1)--------(i)--------(i+1)  --- 
        |           |           |
        |           |           |
---  (i-m-1)------(i-m)------(i-m+1) ---
        |           |           |

the matrix is: 
$$\frac{1}{h^2}[diag(4) + diag(-1,1) + diag(-1,-1) + diag(-1,m) + diag(-1,-m)]$$
simply because 
\begin{align*}-Δu(x,y) &= -∂_{xx}u-∂_{yy}u \\ &=
 \frac{-u(x_{i-1})+2u(x_i)-u(x_{i+1})}{h^2} + \frac{-u(x_{i-m})+2u(x_i)-u(x_{i+m})}{h^2}
\end{align*}

If they understand everything about FDM, you can also introduce FEM to them. That can be done within 2 hours.
Using bilinear Finite Elements on the same lexicographically ordered grid, and the "tensorproduct-trapezoidal-rule" 
$$Q_t(f):=\frac{|T|}{4}\sum_{i=1}^4f(a_i),\qquad a_i \text{ corner points of cell } T,$$ to evaluate the integrals $∫_{T}ϕ^iϕ^jd(x,y)$, will result in the same matrix as 2D FDM (except of the $h^{-2}$-factor).
