Although the duplicate question has an excellent top-voted answer, it may be a bit too dense for beginners. Here's a simple but instructive example. Fibonacci numbers! We all know them, don't we? =)
$\def\nn{\mathbb{N}}$
Let $F_0 = 0$ and $F_1 = 1$ and $F_{n+2} = F_{n+1} + F_n$ for any $n \in \nn$.
$F$ is of course the Fibonacci sequence. To capture the underlying structure of the sequence, we would like to 'factor out' the "$n$" from the above defining recurrence relation. How do we do that? Matrices are a natural solution, as we shall see.
A matrix can be used to encode a linear transformation. How do we view the recurrence relation as a transformation? Our objective is to pass along enough information so that we can generate the sequence by iterating some transformation. Clearly then we need to maintain a pair of consecutive terms, and the transformation is to go from that pair of terms to the next pair:
$(F_{n+1},F_n) \mapsto (F_{n+2},F_{n+1})$.
That immediately gives us the transformation:
$F_{n+2} = 1F_{n+1} + 1F_n$.
$F_{n+1} = 1F_{n+1} + 0F_n$.
Why have I written it like this? Because now we can 'factor out those "$n$"s by fiat! To do so, define:
$\pmatrix{a&b\\c&d} \pmatrix{x\\y} = \pmatrix{ax+by\\cx+dy}$.
Note that the notation on the left is purely notation (for now), and has absolutely nothing to do with matrices or multiplication or anything at all. It just means what is on the right!
Then we can write our transformation as:
$\pmatrix{F_{n+2}\\F_{n+1}} = \pmatrix{1&1\\1&0} \pmatrix{F_{n+1}\\F_{n}}$.
Great! We seem to have 'factored out' the recurrence structure from the terms in the sequence. But as mentioned this is merely notation so far. For this representation to be useful, we need to interpret the notation better.
Firstly, we can treat $\pmatrix{a&b\\c&d}$ as a function on $2$-vectors (which are those things like $\pmatrix{x\\y}$ that consists of a single column of $2$ numbers), specifically:
The function $\pmatrix{a&b\\c&d}$ applied to input $\pmatrix{x\\y}$ gives output $\pmatrix{ax+by\\cx+dy}$.
For example:
$\pmatrix{F_{n+3}\\F_{n+2}} = \pmatrix{1&1\\1&0} \left( \pmatrix{1&1\\1&0} \pmatrix{F_{n+1}\\F_{n}} \right) = \left( \pmatrix{1&1\\1&0} \circ \pmatrix{1&1\\1&0} \right) \pmatrix{F_{n+1}\\F_{n}}$.
Now we can fully capture the iteration of a transformation, which is just repeated application of a function. As usual, for any function $f$ we use "$f^k$" to mean its $k$-fold composition.
The Fibonacci sequence is now completely expressed by:
$\pmatrix{F_{n+1}\\F_n} = \pmatrix{1&1\\1&0}^n \pmatrix{F_{1}\\F_{0}}$.
Notice two important things. Firstly, this representation cleanly divides into the recurrence and the initial values, and works for any choice of $(F_0,F_1)$. Secondly, all we need to do now is to be able to compute $\pmatrix{1&1\\1&0}^n$ directly before applying the resulting function to $\pmatrix{F_{1}\\F_{0}}$. For this we must find out what is the composition of two functions $\pmatrix{a&b\\c&d}$ and $\pmatrix{p&q\\r&s}$. Completely by our above definition we will get:
$\left( \pmatrix{a&b\\c&d} \circ \pmatrix{p&q\\r&s} \right) \pmatrix{x\\y} = \pmatrix{ap+br&aq+bs\\cp+dr&cq+ds} \pmatrix{x\\y}$. [Work this out yourself!]
By the definition of equality of functions we get:
$\left( \pmatrix{a&b\\c&d} \circ \pmatrix{p&q\\r&s} \right) = \pmatrix{ap+br&aq+bs\\cp+dr&cq+ds}$.
Conventionally, each of these three functions is called a matrix, and we shall now define matrix multiplication to be exactly composition! This is the reason for matrix multiplication.
Finally, one should notice that the most efficient way to compute they $n$-fold composition of $\pmatrix{a&b\\c&d}$ is to use the following facts for any function $f$ and $n \in \nn$:
$f^{2n} = (f^n)^2$.
$f^{2n+1} = f^{2n} \circ f$.
This actually yields one of the fastest algorithms possible to compute the $n$-th Fibonacci number with arbitrary precision. (The closed form formula is of course faster for fixed precision.)