Intuition behind PCA eigenvectors
Since the students know about eigenvectors and eigenvalues, I will assume that they know about Lagrangian multipliers as well.

Let's start with computing the variance of data $\vec{x}$ along a given direction $\vec{e}_a$. Since we only care about the direction, $\vec{e}_a$ is a vector of unit length, i.e. 
$$\vec{e}_a^T \vec{e}_a = 1$$

So first, the component of the data along $\vec{e}_a$ is given by
$$ \vec{e}_a^T \vec{x}$$
The variance $\sigma_a^2$ of this component is given by
$$ \sigma_a^2 = E[(\vec{e}_a^T \vec{x} - \vec{e}_a^T \vec{\mu})^2] $$
where $\vec{\mu} = E[\vec{x}]$.
If we re-write this expression a little bit
$$ \sigma_a^2 \\
= E[(\vec{e}_a^T \vec{x} - \vec{e}_a^T \vec{\mu})(\vec{e}_a^T \vec{x} - \vec{e}_a^T \vec{\mu})] \\
= E[\vec{e}_a^T(\vec{x} - \vec{\mu})(\vec{x} - \vec{\mu})^T\vec{e}_a ] \\
= \vec{e}_a^T E[(\vec{x} - \vec{\mu})(\vec{x} - \vec{\mu})^T ] \vec{e}_a \\
= \vec{e}_a^T \Sigma_x \vec{e}_a 
$$
where in the last step we use the covariance matrix $\Sigma_x = E[(\vec{x} - \vec{\mu})(\vec{x} - \vec{\mu})^T ]$.
So the variance along a given direction $\vec{e}_a$ is computed using $\Sigma_x$. So far, so good.

Let's now try to find the direction $\vec{e}_a^*$ where the variance is maximal. Since $\vec{e}_a^*$ is also supposed to only indicate the direction, we demand
$$\vec{e}_a^{*T}\vec{e}_a^{*} = 1 $$
With this we can formulate a constrained optimization problem using the Lagrangian multiplier $\lambda$
$$ L(\vec{e}_a) = \vec{e}_a^T \Sigma_x \vec{e}_a + \lambda (1 - \vec{e}_a^T \vec{e}_a ) $$

Taking the derivative with respect to $\vec{e}_a$ and equating it to zero, we end up with
$$ \frac{\partial L}{\partial \vec{e}_a} = \Sigma_x \vec{e}_a - \lambda \vec{e}_a = 0 $$
i.e.
$$ \Sigma_x \vec{e}_a = \lambda \vec{e}_a $$
This is a cool result. The direction of optimal (maximal) variance is a solution to the eigenvalue problem of $\Sigma_x$. Assuming that we found a normalized eigenvector $\vec{e}_b$, it is easy to interpret the eigenvalue $\lambda$

$$ \vec{e}_b^{T} \Sigma_x \vec{e}_b = \lambda = \sigma_b^{2} $$

So, to summarize:


the variance along a given direction is computed using the covariance matrix
the direction of maximal variance is a normalized eigenvector of the covariance matrix
the eigenvalues are the variances along the direction of the eigenvectors


This means, that the problem of finding the direction of maximal variance reduces to finding the eigenvector corresponding to the maximal eigenvalue. This is about as intuitive as it gets, I believe.

(I skipped the part where we need to show that all eigenvalues are non-negative and thus can be interpreted as variances, but this follows from $\Sigma_x$ being positive semi-definit.)