In general, a b-spline curve will not pass through any of its control points.
There is an example at the bottom of this web page, which explains how repeating knot values will cause a b-spline curve to pass through one of its control points. This technique is typically used with the first and last knots, to force the spline to pass through the first and last control points.
When you repeat knots, this changes the form of the b-spline basis functions, so the ones you cited will not be correct near repeated knots.
To evaluate a cubic b-spline on the interval $[0,1]$, you need a knot sequence that has at least two knot values to the left of 0, and at least two knots to the right of 1. These 6 knots together are needed to define the basis functions that are non-zero on $[0,1]$. So, the knot vector you mentioned $(1,2,3,4, \ldots)$ certainly will not work.
If you want a b-spline curve that passes through a sequence of given points, then you need to use an "interpolating" b-spline. You compute one of these by solving a system of linear equations to get its control points, as explained by @LutzL and on this page.
\begin{align}
f(x) &= a + b (x - x_0) + c (x - x_0)^2 + d (x - x_0)^3
\tag{1}\label{1}
.
\end{align}
Note that \eqref{1} defines $y$ as a function of $x$.
Combined with two values of $x$, $x_0$ and $x_3$,
we have all we need to
make a conversion to equivalent 2D cubic Bezier segment,
defined by its four control points
\begin{align}
P_0(x_0,y_0),\,&P_1(x_1,y_1),\,P_2(x_2,y_2),\,P_3(x_3,y_3)
\tag{2}\label{2}
.
\end{align}
The end points are
\begin{align}
P_0&=(x_0,y_0)=(x_0,f(x_0))=(x_0,a)
\tag{3}\label{3}
,\\
P_3&=(x_3,y_3)=(x_3,f(x_3))
\tag{4}\label{4}
.
\end{align}
2D cubic Bezier segment is defined as usual,
\begin{align}
B_3(t)&=(x(t),y(t))
\tag{5}\label{5}
,\\
x(t)&=x_0\,(1-t)^3+3\,x_1\,(1-t)^2\,t+3\,x_2\,(1-t)\,t^2+x_3\,t^3
\tag{6}\label{6}
,\\
y(t)&=y_0\,(1-t)^3+3\,y_1\,(1-t)^2\,t+3\,y_2\,(1-t)\,t^2+y_3\,t^3
,\quad t\in[0,1]
\tag{7}\label{7}
.
\end{align}
We know that $x$ is linear in $t$, so
we must have
\begin{align}
x_0\,(1-t)+x_3\,t
&=x_0\,(1-t)^3+3\,x_1\,(1-t)^2\,t+3\,x_2\,(1-t)\,t^2+x_3\,t^3
\tag{8}\label{8}
,\\
(x_3-x_0)t+x_0
&=(3x_1-x_0-3x_2+x_3)t^3+(3x_0+3x_2-6x_1)t^2+3(x_1-x_0)t+x_0
\quad \forall t\in[0,1]
\tag{9}\label{9}
,
\end{align}
hence we have $x_1,x_2$ evenly distributed between the endpoints $x_0$ and $x_3$:
\begin{align}
x_1 &= \tfrac13 (2x_0+x_3)=x(t)\Big|_{t=1/3}
\tag{10}\label{10}
,\\
x_2 &= \tfrac13 (x_0+2x_3)=x(t)\Big|_{t=2/3}
\tag{11}\label{11}
.
\end{align}
Corresponding $y$ points on the curve \eqref{1}
are
\begin{align}
y(\tfrac13)&=f(x_1)
\tag{12}\label{12}
,\\
y(\tfrac23)&=f(x_2)
\tag{13}\label{13}
.
\end{align}
The last pair of equations is a linear system
with two unknowns, $y_1$ and $y_2$
which can be trivially solved as
\begin{align}
y_1
&=
\tfrac13\,b\,(x_3-x_0)+a
\tag{14}\label{14}
,\\
y_2
&=
\tfrac13(x_3-x_0)(c(x_3-x_0)+2b)+a
\tag{15}\label{15}
.
\end{align}
Example
\begin{align}
a &= 7
,\quad
b = 2
,\quad
c = 2
,\quad
d = -1
\tag{16}\label{16}
,\\
x_0 &= -1
,\quad
x_3 = 3
\tag{17}\label{17}
,\\
y_0&=a=7
,\quad
y_3=f(3)=-17
\tag{18}\label{18}
,\\
x_1 &= \tfrac13(2x_0+x_3)=\tfrac13
\tag{19}\label{19}
,\\
x_2 &= \tfrac13(x_0+2x_3)=\tfrac53
\tag{20}\label{20}
,\\
y_1
&=
\tfrac13\,b\,(x_3-x_0)+a
=\tfrac{29}3
\tag{21}\label{21}
,\\
y_2
&=
\tfrac13(x_3-x_0)(c(x_3-x_0)+2b)+a
=23
\tag{22}\label{22}
.
\end{align}
Best Answer
If you want to use the points $\mathbf{P}_0, \ldots, \mathbf{P}_n$ as control points of a spline curve, then the curve is just $$ \mathbf{C}(t) = \sum_{i=0}^n B_i(t)\mathbf{P}_i $$ where the $B_i$ are b-spline basis functions. This formula allows you to compute points on the curve, which you can then use to draw it. Note that, in general, this curve will not pass through the points $\mathbf{P}_0, \ldots, \mathbf{P}_n$.
If you want your curve to pass through the points $\mathbf{P}_0, \ldots, \mathbf{P}_n$, then you first need to construct an interpolating spline.
You will need to use a parametric spline that gives $x$ and $y$ as functions of some parameter $t$. To construct the curve, you must first assign a parameter value $t_i$ to each point $P_i$. These parameter values are entirely fabricated -- they are not an intrinsic part of the original interpolation problem. Furthermore, different choices for these parameter values will produce significantly different curves. The usual approach is to choose them so that $t_{i+1}-t_i$ is equal to the chord-length between the points $\mathbf{P}_i$ and $\mathbf{P}_{i+1}$.
You solve systems of linear equations to get the $x$ and $y$ coordinates of control points $\mathbf{Q}_0, \ldots, \mathbf{Q}_n$. Then the equation of the curve is again $$ \mathbf{C}(t) = \sum_{i=0}^n B_i(t)\mathbf{Q}_i $$
Note that we haven't really made use of the fact that the points $\mathbf{P}_i$ or $\mathbf{Q}_i$ are 2D; this same construction will work in 3D, or in fact in any finite-dimensional space.
A good reference is this web site.