We can also use Linear Algebra to solve
Ordinary Differential Equations
An ODE is of the form
$$\underline{u}'(t) = A \underline{u}(t) + \underline{b}(t)$$
with $A \in \mathbb{C}^{n \times n}$ and $\underline{b}(t) \in \mathbb{C}^{n \times 1}$. If we have an initial condition
$$\underline{u}(t_0) = \underline{u_0}$$
this is an initial value problem. Assuming the entries of $\underline{b}(t)$ are continuous on $[t_0,T]$ for some $T > t_0$, Picard-Lindelöf provides a unique solution on that interval. If $A$ is diagonalisable, the solution of the homogeneous initial value problem is easy to compute.
Let
$$P^{-1} A P = \Lambda = \operatorname{diag}(\lambda_1, \dots, \lambda_n),$$
where $P = \begin{pmatrix} x_1 & \dots & x_n \end{pmatrix}$. Defining $\tilde{\underline{u}}:= P^{-1} \underline{u}(t)$ and $\tilde{\underline{u_0}} = P^{-1} \underline{u_0}$, the IVP reads
$$\tilde{\underline{u}}'(t) = \Lambda \tilde{\underline{u}}(t), \; \tilde{\underline{u}}(t_0) = \tilde{\underline{u_0}} =: \begin{pmatrix} c_1 & \dots & c_n \end{pmatrix}^T.$$
These are simply $n$ ordinary, linear differential equations
$$\tilde{u_j}'(t) = \lambda_j \tilde{u_j}(t), \; \tilde{u_j}(t_0) = c_j$$
for $j = 1, \dots, n$ with solutions $\tilde{u_j}(t) = c_j e^{\lambda_j(t-t_0)}$. We eventually retrieve $\underline{u}(t) = P \tilde{\underline{u}}(t)$.
Example: We can write
$$x''(t) = -\omega^2 x(t), \; x(0) = x_0, \; x'(0) = v_0$$
as $\underline{u}'(t) = A \underline{u}(t), \; \underline{u}(0) = \underline{u_0}$, where $\underline{u}(t) = \begin{pmatrix} x(t)&x'(t) \end{pmatrix}^T$ and
$$A = \begin{pmatrix} 0&1\\ -\omega^2&0 \end{pmatrix} \text{ and } \underline{u_0} = \begin{pmatrix} x_0\\ v_0 \end{pmatrix}.$$
Computing eigenvalues and eigenvectors, we find
$$\underline{u}(t) = c_1 e^{i \omega t} \begin{pmatrix} 1\\ i \omega \end{pmatrix} + c_2 e^{-i \omega t} \begin{pmatrix} 1 \\ -i \omega \end{pmatrix}.$$
Using the initial condition, we find $x(t) = x_0 \cos(\omega t) + \frac{v_0}{\omega} \sin(\omega t)$.
Matrix exponential: I don't know if your students already are familiar with the matrix exponential, but using it we find a solution of the homogeneous initial value problem to be given by
$$\underline{u}(t) = e^{(t-t_0)A} \underline{u_0}.$$
To solve the inhomogeneous differential equation, we use can vary the constants. Since every solution of the homogeneous system $\underline{u}'(t) = A \underline{u}(t)$ is of the form $\underline{u}(t) = e^{tA} \underline{c}$ for some constant vector $\underline{c}$, we set $\underline{u_p}(t) = e^{tA} \underline{c}(t)$ and find by plugging in
$$\underline{c}'(t) = e^{-tA} \underline{b}(t).$$
Thus,
$$\underline{u}(t) = e^{(t-t_0)A} \underline{u_0} + \int_{t_0}^t e^{(t-s)A} \underline{b}(s) \, \mathrm ds.$$