The undergrad course I took which included B-splines spent a lot of time first on Bézier curves. You might not necessarily want to spend much time on them, but I think they can motivate a variant on method 1 by defining Bézier curves with de Casteljau's algorithm and B-spline curves with de Boor's algorithm. This assumes that you're at least as interested in the curves as in the basis functions, but since you talk about applications...
Note that I haven't tested this on anyone (undergrad or no) and the latter half has somewhat the flavour of experimental mathematics, which may not suit your pedagogical style.
Basics
We start from the concept of a spline curve as a curve defined by some control points and weighting functions, and we parameterise the weighting functions over an interval of the real line. So generically we have $P(t) = \sum_k w_k(t) P_t$ where $w_k$ is the $k$th weighting function and $P_k$ is the $k$th control point. One highly desirable property of the weighting functions is that for all $t$ in the interval on which the line is defined, $\sum_k w_k(t) = 1$, because this is necessary and sufficient for an affine transformation of the control points to transform the entire line consistently.
Recursive interpolation
A straightforward way to achieve this unit weight property is to choose weighting functions corresponding to recursive interpolation. In other words, the evaluation process is to define a series of sequences of control points starting with the original control points ($P_{0,k} = P_k$) and repeatedly interpolating between adjacent control points, reducing the number of points by one each time, until we get down to a single point. In general, $$P_{i,k} = (1 - f_{i,k}(t)) P_{i-1,k} + f_{i,k}(t) P_{i-1,k+1}$$
(This may be a good point at which to show that if the $f_{i,k}$ all map into $[0,1]$ then we get the convex hull property).
Example 1: Bézier
The simplest such scheme is to take $[0,1]$ as our interval of the real line and $f_{i,k}(t) = t$. Then by induction we can show that the corresponding weighting functions are $w_k(t) = \binom{n}{k} (1-t)^{n-k} t^k$ when we have $n+1$ control points numbered $P_0$ to $P_n$. This is the Bézier spline, which is very popular but has some important drawbacks. In particular, we can see that each weighting function has support over the entirety of $(0,1)$. We can attempt to fix this and achieve local control by chaining a series of Bézier curves, which we now treat as segments in a larger curve. For simple continuity we require the last control point of one segment to be the first control point of the next segment; if we want higher level continuity then we impose conditions on more control points. (There's an elegant umbral calculus approach to Bézier curves as $P(t) = (1 - t + t \mathcal{P})^n$ where $\mathcal{P}^k$ corresponds to $P_k$ which allows us to derive the condition for $c$th order continuity from segment $P$ to subsequent segment $Q$ as $$\forall j \le c: Q_j = 2^j \sum_{k=0}^j (-2)^k \binom{j}{k} P_{n-k}$$
so that if we want a degree of continuity at the joins which is half of the degree of the segments then moving a control point in one segment can have knock-on effects beyond the adjacent segment; but I suspect that this level of detail is too much of a digression for you. It's more straightforward to show that the $c$th derivative of segment $P$ at $P_n$ depends on the last $c+1$ control points in the segment, and similarly the $c$th derivative of segment $Q$ at $Q_0$ depends on the first $c+1$ control points, so that equating them necessarily establishes constraints between those control points).
Example 2: cardinal B-spline
A natural solution is to say that if we're going to have constraints between the control points which directly affect adjacent segments of our curve, we should take the simplest possible constraints: identification. Therefore we'll have a single set of control points $P_k$ and we'll say that for $t \in [a, a+1)$ we'll perform a recursive interpolation involving points $P_a$ to $P_{a+n}$. (This may not be the standard indexing, for which I apologise). Still in the interests of simplicity, we'll take the $f_{i,k}$ as linear functions: $f_{i,k}(t) = \alpha_{i,k} t + \beta_{i,k}$. And in the interests of uniformity we'll refine the interpolation to take into account that we're handling multiple segments, and update it to $$P_{i,k} = (1 - f_{i,k}(t-a)) P_{i-1,k} + f_{i,k}(t-a) P_{i-1,k+1}$$
Then the question is what values of $\alpha_{i,k}$ and $\beta_{i,k}$ maximise the continuity between the segments without imposing constraints between the $P_k$.
In $[a, a+1)$ we have (introducing some hopefully transparent notation for brevity)
- $P_{0,k} = P_{a+k}$ for $0 \le k \le n$.
- $P_{1,k} = \overline{f_{1,k}} P_{a+k} + f_{1,k} P_{a+k+1}$ for $0 \le k \le n-1$.
- $P_{2,k} = \overline{f_{1,k}} \overline{f_{2,k}} P_{a+k} + (f_{1,k} \overline{f_{2,k}} + \overline{f_{1,k+1}} f_{2,k}) P_{a+k+1} + f_{1,k+1} f_{2,k} P_{a+k+2}$ for $0 \le k \le n-2$
etc.
In the interests of having a digestible example, consider the case $n=2$.
$P(t) = P_{n,0} = \overline{f_{1,0}} \overline{f_{2,0}} P_{a} + (f_{1,0} \overline{f_{2,0}} + \overline{f_{1,1}} f_{2,0}) P_{a+1} + f_{1,1} f_{2,0} P_{a+2}$
In the segment $[a, a+1)$,
$$\lim_{t \to a+1} P(t) = \overline{f_{1,0}} \overline{f_{2,0}} (1) P_{a} + (f_{1,0} \overline{f_{2,0}} + \overline{f_{1,1}} f_{2,0}) (1) P_{a+1} + f_{1,1} f_{2,0} (1) P_{a+2}$$ and in the segment $[a+1, a+2)$ we have
$$P(a+1) = \overline{f_{1,0}} \overline{f_{2,0}} (0) P_{a+1} + (f_{1,0} \overline{f_{2,0}} + \overline{f_{1,1}} f_{2,0}) (0) P_{a+2} + f_{1,1} f_{2,0} (0) P_{a+3}$$
So for simple continuity without imposing constraints on the control points we require
$$\begin{eqnarray*}
(1 - \alpha_{1,0} - \beta_{1,0})(1 - \alpha_{2,0} - \beta_{2,0}) &=& 0 \\
(\alpha_{1,0} + \beta_{1,0}) (1 - \alpha_{2,0} - \beta_{2,0}) + (1 - \alpha_{1,1} - \beta_{1,1}) (\alpha_{2,0} + \beta_{2,0}) &=& (1 - \beta_{1,0})(1 - \beta_{2,0}) \\
\alpha_{1,1} \alpha_{2,0} + \alpha_{1,1} \beta_{2,0} + \alpha_{2,0} \beta_{1,1} &=& \beta_{1,0} - \beta_{1,0} \beta_{2,0} + \beta_{2,0} \\
0 &=& \beta_{1,1} \beta_{2,0}
\end{eqnarray*}$$
This still leaves some degrees of freedom, so we can consider continuity of the first derivative, which (some algebra-bashing omitted) yields
\begin{eqnarray*}
\alpha_{1,0}(1 - 2 \alpha_{2,0} - \beta_{2,0}) + \alpha_{2,0} (1 - \beta_{1,0}) &=& 0 \\
\alpha_{1,0} (2 - 2 \alpha_{2,0} - 2 \beta_{2,0}) - \alpha_{1,1} (2 \alpha_{2,0} + \beta_{2,0}) + \alpha_{2,0} (2 - 2 \beta_{1,0} - \beta_{1,1}) &=& 0 \\
\alpha_{1,0} (1 - \beta_{2,0}) - \alpha_{1,1} (2 \alpha_{2,0} + 2 \beta_{2,0}) + \alpha_{2,0} (1 - \beta_{1,0} - 2 \beta_{1,1}) &=& 0 \\
\alpha_{1,1} \beta_{2,0} + \alpha_{2,0} \beta_{1,1} &=& 0
\end{eqnarray*}
and combining the two sets of continuity constraints we get $\alpha_{1,0} = \tfrac12$, $\alpha_{1,1} = \tfrac12$, $\alpha_{2,0} = 1$, $\beta_{1,0} = \tfrac12$, $\beta_{1,1} = 0$, $\beta_{2,0} = 0$.
If we now suggest considering $\alpha_{i,k} = \frac1{n+1-i}$, $\beta_{i,k} = \frac{n-k-i}{n+1-i}$ then it's no longer entirely unmotivated...
The full B-spline
Thus far we haven't mentioned knots, but there's no fundamental reason why we should segment the parameter space at the integers...
Once the full generality has been introduced, and since we came via Bézier curves, it's probably worth throwing out a comment about the Bézier basis functions corresponding to B-spline basis functions with knot vectors $0^u 1^v$.
