I am aware of two main ways for deriving the intersection between a ray and a triangle. I believe most ray-tri intersection routines are variants of those. In fact Möller–Trumbore is not something special from a purely mathematical perspective, it's just a specific grouping of terms when solving a 3x3 system.
The first derivation involves finding barycentric and ray coordinates simultaneously which results in a 3x3 system. The second chooses to first intersect the ray with the plane, in which the triangle lies, and then finds the barycentric coordinates.
0.1) Parametric Ray Definition
A parametric definition of a line given an origin $\vec{o}$ and a direction $\vec{d}$ is:
$$\vec{r}(t) = \vec{o} + t\vec{d}.$$
A ray is the set of points on one side of the origin of the line:
$$\mathcal{R} = \{\vec{r}(t) = \vec{o} + t\vec{d}\,:\,t\in[0,\infty)\}.$$
You should convince your students that this is indeed a ray. A good approach would be to draw several point on the ray for different $t$ vaues (e.g. 0, 1, 0.5, but also note that -1 is on the line but not part of the ray).
0.2) Parametric Triangle Definition
A parametric definition of a plane given a point on it $\vec{v}_0$ and two vectors $\vec{f}_1 = \vec{v}_1 - \vec{v}_0$, $\vec{f}_2 = \vec{v}_2-\vec{v}_0$ is:
$$\vec{v}(\beta, \gamma) = \vec{v}_0 + \beta\vec{f}_1 + \gamma\vec{f}_2$$
The pair $(\beta, \gamma)$ are the coordinates of $\vec{v}$ in the 2D coordinate system with origin $\vec{v}_0$ and basis vectors $\vec{f}_1, \vec{f}_2$. For those to be basis vectors we additionally require that $\vec{f}_1, \vec{f}_2$ are not parallel/linearly dependent (i.e. the triangle is non-degenerate).
You should convince the students that $\vec{v}(\beta, \gamma)$ indeed spans the whole plane for different $(\beta, \gamma)$. For instance you could draw the points $(0,0), (1,0), (0,1), (1,1)$ which are the vertices of a parallelogram. Additionally, the points forming the 3 edges of the triangle are: $(0,\gamma); \gamma \in[0,1]$, $(\beta, 0); \beta \in [0,1]$, and $(\beta, \gamma); \beta \geq 0, \gamma \geq 0, \beta + \gamma = 1$. Then the set of points forming the triangle are:
$$\mathcal{T} = \{\vec{v}(\beta, \gamma)\,:\,0\leq \beta, 0\leq \gamma, \beta+\gamma \leq 1\}.$$
This should definitely be accompanied with a drawing of the geometry so that the students may grasp its meaning.
0.2.1) Barycentric Coorindates
This is something that I would personally leave for the very end, since it is not integral to any of the derivation, is less general, and more complex.
The coordinates $\beta, \gamma$ are in fact the barycentric coordinates of $\vec{v}$ associated with $\vec{v}_1, \vec{v}_2$. In fact we can rewrite our formulation as:
$$\begin{align}\vec{v}(\beta, \gamma) &= \vec{v}_0 + \beta \vec{f}_1 + \gamma \vec{f}_2 \\ \vec{v}(\beta, \gamma) &= \alpha\vec{v}_0 + \beta\vec{v}_1 + \gamma\vec{v}_2, \alpha = 1-\beta -\gamma,\end{align}$$
where $\alpha, \beta, \gamma$ are the barycentric coordinates of point $\vec{v}$ with respect to the triangle $\vec{v}_0, \vec{v}_1, \vec{v}_2$. The choice of $(\beta, \gamma)$ over $(\alpha, \beta)$ or $(\alpha, \gamma)$ was arbitrary. Indeed using $\alpha+\beta+\gamma=1$ I could have equivalently written:
$$\vec{v}(\alpha, \beta) = \vec{v}_2 + \alpha (\vec{v}_0-\vec{v}_2) + \beta (\vec{v}_1-\vec{v}_2) \\ \text{or} \\
\vec{v}(\alpha, \gamma) = \vec{v}_1 + \alpha (\vec{v}_0-\vec{v}_1) + \gamma (\vec{v}_2-\vec{v}_1).
$$
However, the fact that $\alpha, \beta, \gamma$ are the barycentric coordinates is not really important, and can be kept as an additional fact in order to not overload the students. The 2D basis interpretation is simpler and much more general.
0.3) Problem Formulation
We want to find the closest intersection (the one with smallest $t$) between $\mathcal{R}$ and $\mathcal{T}$. Let $\vec{p}$ be an intersection point, i.e. $\vec{p} \in \mathcal{R}$ and $\vec{p} \in \mathcal{T}$, then there exist $t\geq 0, \beta\geq 0, \gamma\geq 0, \beta+\gamma\leq 1$ such that:
$$\vec{o}+t\vec{d} = \vec{p} = \vec{v}_0 + \beta\vec{f}_1 + \gamma\vec{f}_2.$$
Conversely, if those do not exist then we'll say that there is no intersection. Beyond this point the only difference is how one groups the terms when solving the above $3\times 3$ system. For example Möller–Trumbore is just a special case of a specific solution order for the unknowns $t, \beta, \gamma$.
1) Simultaneous Solution
Let us rewrite the equations in matrix form:
$$\vec{o}+t\vec{d} = \vec{v}_0 + \beta\vec{f}_1 + \gamma\vec{f}_2$$
$$\begin{bmatrix}\vec{d} & \vec{f}_1 & \vec{f}_2 \end{bmatrix} \begin{bmatrix} -t \\ \beta \\ \gamma \end{bmatrix} = \vec{o} - \vec{v}_0.$$
It should be made clear for the students that $M = \begin{bmatrix}\vec{d} & \vec{f}_1 & \vec{f}_2 \end{bmatrix}$ is a 3x3 matrix, didactically, being explicit there may help:
$$M = \begin{bmatrix}d_x & f_{1,x} & f_{2,x} \\ d_y & f_{1,y} & f_{2,y} \\ d_z & f_{1,z} & f_{2,z}\end{bmatrix}.$$
Solving the system is as simple as:
$$\begin{bmatrix} -t \\ \beta \\ \gamma \end{bmatrix} = M^{-1}(\vec{o}-\vec{v}_0)$$
Either $3\times 3$ matrix inversion followed by matrix-vector multiplication can be used, or the equivalent Cramer's rule formulation. There's also a trick involving the cross product that works for $3\times 3$ matrix inversion that some people like to use (which is elegant, but didactically not very useful):
$$M^{-1} = \frac{1}{det[M]}adj(M) = \frac{1}{\vec{d}\cdot(\vec{f}_1\times\vec{f}_2)}\begin{bmatrix} (\vec{f}_1\times\vec{f}_2)^T \\ (\vec{f}_2 \times \vec{d})^T \\ (\vec{d}\times \vec{f}_1)^T\end{bmatrix}$$
You can identify the terms in Möller–Trumbore's algorithm with the terms from above, albeit I would altogether just recommend skipping Möller–Trumbore, and giving your students a 3x3 matrix inverse procedure and a matrix-vector multiplication procedure which should be enough to find the vector $\begin{bmatrix} -t & \beta & \gamma \end{bmatrix}^T$.
Once $t, \beta, \gamma$ are computed one needs to check that this is indeed a valid intersection (i.e. it belongs to both $\mathcal{R}$ and $\mathcal{T}$):
$$(0\leq t < \infty) \land (0\leq \beta) \land (0\leq \gamma) \land (\beta+\gamma\leq 1).$$
In C++ the above check would fail also when the matrix is not invertible, since NaNs and infinities will fail this check, so no additional code is necessary for the special cases. It is instructive however to consider what the special cases represent. One option is that $\vec{f}_1$ and $\vec{f}_2$ are parallel/linearly dependent - that means that we have a degenerate triangle (it's rather a line segment in this case). Another option is that $\vec{d}, \vec{f}_1, \vec{f}_2$ are coplanar/linearly dependent. This may occur if $\vec{d}$ is parallel to the plane spanned by $\vec{f}_1,\vec{f}_2$. Then if $\vec{o}$ doesn't lie in the plane there is no intersection. If $\vec{o}$ lies in the plane, then there are infinitely many intersections, and the closest one is at $t=0$. For convenience one can treat this as a non-intersection however, since the ray is inside the plane, which is infinitely thin.
2) Split Solution
The other solution method that I like and which has a nice geometric interpretation is the one where we first intersect the plane in which the triangle lies, and only after do we find the $\beta, \gamma$ coordinates. This is loosely equivalent to solving for the $t$ variable in the system, and then substituting back in order to get a $2\times 2$ system. We know that the plane has a normal $\vec{n} = \vec{f}_1\times\vec{f}_2$ and that the set of points on the plane are:
$$\mathcal{P} = \{\vec{p}\,:\, (\vec{p}-\vec{v}_0)\cdot \vec{n} = 0\}.$$
Then plugging in $\vec{r}(t)$ for $\vec{p}$ yields:
$$(\vec{o}+t\vec{d}-\vec{v}_0)\cdot \vec{n} = 0 \implies t = \frac{(\vec{v}_0-\vec{o})\cdot\vec{n}}{\vec{d}\cdot \vec{n}}.$$
Note also that this is exactly what you'll get using Cramer's rule for $t$ using the $3\times 3$ system. Provided that $t \in [0,\infty)$ we can use the fact that we know the intersection point on the plane: $\vec{p} = \vec{o}+t\vec{d}$. We are now left with 3 equations but two unknowns:
$$\vec{p} = \vec{v}_0 + \beta\vec{f}_1 + \gamma\vec{f}_2.$$
Ideally we would like to reduce this to a system of $2$ equations, this can be done by dotting the equations with $\vec{f}_1$ and $\vec{f}_2$:
$$\vec{f}_1\cdot(\vec{p}-\vec{v}_0) = \beta\vec{f}_1\cdot\vec{f}_1 + \gamma\vec{f}_1\cdot\vec{f}_2 \\
\vec{f}_2\cdot(\vec{p}-\vec{v}_0) = \beta\vec{f}_2\cdot\vec{f}_1 + \gamma\vec{f}_2\cdot\vec{f}_2.$$
It's also instructive to look at what this looks like in terms of matrix notation:
$$F = \begin{bmatrix} \vec{f}_1 & \vec{f}_2\end{bmatrix} \\
\vec{p}-\vec{v}_0 = F \begin{bmatrix} \beta \\ \gamma \end{bmatrix} \\
F^T(\vec{p}-\vec{v}_0) = F^TF\begin{bmatrix} \beta \\ \gamma \end{bmatrix}.$$
Then the solution reads:
$$\begin{bmatrix} \beta \\ \gamma \end{bmatrix} = (F^TF)^{-1}F^T(\vec{p}-\vec{v}_0),$$
and one must check that $(0\leq \beta) \land (0\leq \gamma) \land (\beta+\gamma\leq 1)$ as before.
2.1) Algebraic and Geometric Meaning
These are additional mathematically interesting details, that are probably best left for the very end.
The matrix $F^TF$ is a $2\times 2$ matrix that is trivial to invert and has the following form:
$$F^TF = \begin{bmatrix} \vec{f}_1\cdot\vec{f}_1 & \vec{f}_1\cdot\vec{f}_2 \\ \vec{f}_1\cdot\vec{f}_2 & \vec{f}_2\cdot\vec{f}_2\end{bmatrix}.$$
The matrix $G=F^TF$ is in fact the metric tensor on the plane. To motivate the metric tensor you can describe it as a tool to measure angles and lengths in a non-orthonormal coordinate system. Given two vectors $\vec{u}$ and $\vec{w}$ in the plane described through $\vec{u}_P = \begin{bmatrix} \beta_u & \gamma_u\end{bmatrix}^T$ and $\vec{w}_P = \begin{bmatrix} \beta_w & \gamma_w\end{bmatrix}^T$ the metric tensor can be used to compute the equivalent of the dot product between those:
$$\vec{u}_P^TG\vec{w}_P = \vec{w}_P^TG\vec{u}_P = \vec{u}\cdot\vec{w} = \|\vec{u}\|\|\vec{w}\|\cos\angle(\vec{u},\vec{w}).$$
In the above $\vec{u}, \vec{w}$ are the representations of the vectors with respect to the canonical basis $(1,0,0),\, (0,1,0),\, (0,0,1)$. It is also clear that for an orthonormal basis $\vec{f}_1, \vec{f}_2$, the metric tensor is the identity matrix, recovering the usual dot product. In fact any inner product in a finite-dimensional real space has the form: $\vec{u}^TG\vec{w}$ for $G$ symmetric and positive definite. So this can be used to introduce both the concept of metric tensor and inner product with a suitable motivation. Using the definition for length: $\|\vec{u}_P\|^2 = \vec{u}_P^TG\vec{u}_P = \vec{u}^T\vec{u}$ it is also clear that the metric tensor helps with measuring length in non-orthonormal coordinate systems.
Another interpretation of the system is that:
$$F^T(\vec{p}-\vec{v}_0) = F^TF\begin{bmatrix} \beta \\ \gamma \end{bmatrix}$$
are the normal equations, which are used when we have an overdetermined system (more linearly independent equations than unknowns, note that $F$ is rectangular with more rows than columns). In fact even if $\vec{p}$ was not in the plane of the triangle, the above would find the $\beta, \gamma$ for its projection on this plane, since the solution of the normal equations is a minimizer of:
$$\left\|\vec{v}_0 + F\begin{bmatrix} \beta \\ \gamma \end{bmatrix} - \vec{p}\right\|^2.$$
Moreover in:
$$\begin{bmatrix} \beta \\ \gamma \end{bmatrix} = (F^TF)^{-1}F^T(\vec{p}-\vec{v}_0),$$
the matrix $(F^TF)^{-1}F^T$ is known as the Moore–Penrose (left pseudo-)inverse. So there are in fact many algebraic and geometric interpretations that can be made in order to connect this problem to topics that the students have studied or are about to study. It's a good motivational tool.
