The reason your answer differs is because of spline boundary conditions. Consider the spline $B_{j,1}$. To compute this we need $B_{j,0}$, $B_{j+1,0}$, and $\omega_{j,1}$ and $\omega_{j+1,1}$. $$ B_{j,0} = \begin{cases}1 & x \in [t_j,t_{j+1}]\\0&\mathrm{else}\end{cases}\\ B_{j+1,0} = \begin{cases}1 & x \in [t_{j+1},t_{j+2}]\\0&\mathrm{else}\end{cases}\\ \omega_{j,1} = \frac{x - t_{j}}{t_{j+1}-t_j}\\ \omega_{j+1,1} = \frac{x - t_{j+1}}{t_{j+2}-t_{j+1}} $$ Notice that to compute $B_{j,1}$ in the interval $[t_{i},t_{i+1}]$ requires knowing the knots $t_{i}$ through $t_{i+2}$. If we try evaluating in $[t_{4},t_{5}]$ we need to know $t_6$, which is undefined. Additionally, there are a few weird things which happen when we try to evaluate in $[t_0,t_1]$. There could be some $t_{-1}$ which would affect the result in this interval. In general to compute $B_{j,k}$ we only know for certainty what its value is between $[t_{k},t_{n-1-k}]$.
We can adjust our original knots vector by adding extra knots $t'$ on both sides to allow us to evaluate our spline for the full original interval $[t_0,t_{n-1}]$. The most common way is to simply duplicate the end knots $k$ times on each side, so we have $t'=[t_0,\ldots,t_0, \ldots,t_{n-1},\ldots,t_{n-1}]$. This is what the page you linked does. Some other options are what the extrapolate option does in SciPy: for extrapolate='periodic' the knots vector is considered to be periodic. For extrapolate=True the resultant polynomial is simply extrapolated as-is outside the interval $[t_{k},t_{n-1-k}]$.
Here's code for duplicating the figures in the linked page:
import matplotlib.pyplot as plt import numpy as np from scipy.interpolate import BSpline x = np.linspace(-1,10,1000) t = np.array([1, 2, 3, 4, 5, 6]) k = 1 tp = np.concatenate((np.repeat(t[0], k), t, np.repeat(t[-1], k))) e = np.eye(tp.size) plt.figure() for i in range(t.size-1-k): c = e[i+k] spl = BSpline(tp, c, k, extrapolate=False) y = spl(x) plt.plot(x,y) ```