Substituting variables $(\vec{r}_i,\vec{r}_j,\vec{r}_k) = (A,B,C)$, with all z-coordinates equal to zero, into the 'in-class' formula gives the following:
$$\vec{r}_c = A+\frac{({\lVert C-A \rVert}^2(B-A)-{\lVert B-A \rVert}^2(C-A)) \times (C-A) \times (B-A)}{2{\lVert (C-A) \times (B-A) \rVert}^2}$$
$$\vec{r}_c = A+\frac{(({\lVert A \rVert}^2+{\lVert C \rVert}^2-2A \cdot C)(B-A)-({\lVert A \rVert}^2+{\lVert B \rVert}^2-2A \cdot B)(C-A)) \times (C-A) \times (B-A)}{2{\lVert (C-A) \times (B-A) \rVert}^2} $$
$$\vec{r}_c = A+\frac{\left({\lVert A \rVert}^2 (B-C) + {\lVert B \rVert}^2 (A-C) + {\lVert C \rVert}^2 (B-A) -2\left[(A\cdot C)(B-A)-(A\cdot B)(C-A)\right]\right) \times (C-A) \times (B-A)}{2{\lVert (C-A) \times (B-A) \rVert}^2} $$
Taking a hint from the denominator, I assume the cross product follows the following order. (The terms to the right combine to form a unit vector.)
$$\vec{r}_c = A+\frac{{\lVert A \rVert}^2 (B-C) + {\lVert B \rVert}^2 (A-C) + {\lVert C \rVert}^2 (B-A) -2\left[(A\cdot C)(B-A)-(A\cdot B)(C-A)\right]}{2{\lVert (C-A)\times(B-A) \rVert}} \times \frac{(C-A)\times(B-A)}{\lVert (C-A)\times(B-A)\rVert}$$
$$\vec{r}_c = A+\left(\frac{{\lVert A \rVert}^2 (B-C) + {\lVert B \rVert}^2 (A-C) + {\lVert C \rVert}^2 (B-A)}{2(A_x(B_y-C_y)+B_x(C_y-A_y)+C_x(A_y-B_y))}-\frac{(A\cdot C)(B-A)-(A\cdot B)(C-A)}{\lVert (C-A)\times(B-A)\rVert}\right)\times \hat z$$
Let $Y=B-A$, $Z=C-A$, $A\cdot B=X\cdot Y$, and $A\cdot C=X\cdot Z$, giving the following:
$$\vec{r}_c = A+\left(\frac{{\lVert A \rVert}^2 (B-C) + {\lVert B \rVert}^2 (A-C) + {\lVert C \rVert}^2 (B-A)}{2(A_x(B_y-C_y)+B_x(C_y-A_y)+C_x(A_y-B_y))}-\frac{(X\cdot Z)Y-(X\cdot Y)Z}{\lVert (C-A)\times(B-A)\rVert}\right)\times \hat z$$
$$\vec{r}_c = A+\left(\frac{{\lVert A \rVert}^2 (B-C) + {\lVert B \rVert}^2 (A-C) + {\lVert C \rVert}^2 (B-A)}{2(A_x(B_y-C_y)+B_x(C_y-A_y)+C_x(A_y-B_y))}-\frac{X\times(Y\times Z)}{\lVert (C-A)\times(B-A)\rVert}\right)\times \hat z$$
$$\vec{r}_c = A+\left(\frac{{\lVert A \rVert}^2 (B-C) + {\lVert B \rVert}^2 (A-C) + {\lVert C \rVert}^2 (B-A)}{2(A_x(B_y-C_y)+B_x(C_y-A_y)+C_x(A_y-B_y))}+X\times\frac{(C-A)\times(B-A)}{\lVert (C-A)\times(B-A)\rVert}\right)\times \hat z$$
$$\vec{r}_c = A+(X\times\hat z)\times\hat z +\frac{{\lVert A \rVert}^2 (C-B) + {\lVert B \rVert}^2 (A-C) + {\lVert C \rVert}^2 (B-A)}{2(A_x(B_y-C_y)+B_x(C_y-A_y)+C_x(A_y-B_y))}\times\hat z$$
Now, if I knew how to solve for $-X$, I'd show the following:
$$\vec{r}_c = \frac{{\lVert A \rVert}^2 (C-B) + {\lVert B \rVert}^2 (A-C) + {\lVert C \rVert}^2 (B-A)}{2(A_x(B_y-C_y)+B_x(C_y-A_y)+C_x(A_y-B_y))}\times\hat z$$
In other words, the two formulas are off by a factor of -1, and it appears to be the Wikipedia one that maps improperly, as I had been 'fixing' it by setting yc=-yc -- the x-axis is symmetric, so I didn't notice that it was also flipped.
The 'in-class' formula needs to be implemented like this: (@Chris Taylor hit it on the head!)
C = [I,0]+(cross((norm(K-I))^2*[J-I,0]-(norm(J-I))^2*[K-I,0],cross([K-I,0],[J-I,0])))/(2*(norm(cross([K-I,0],[J-I,0])))^2);
C = C(1,1:2);
I'll have to try it out later, as it's time to head to class!
Here is a step by step approach.
$1.$ As we have the coordinates of $B$ and $C$, we can compute a unit normal vector to the line segment $BC$. This is how yo may do this
$$\begin{align}
{\bf{L}}&= {\bf{x}}_B-{\bf{x}}_C\\
&=(x_B-x_C){\bf{i}}+(y_B-y_C){\bf{j}} \\
\\
{\bf{N}}&= {\bf{k}} \times {\bf{L}}\\
&=-(y_B-y_C){\bf{i}}+(x_B-x_C){\bf{j}} \\
\\
{\bf{N}} \cdot {\bf{L}} &= 0 \\
\\
{\bf{n}} &= \frac{\bf{N}}{\|\bf{N}\|} \\
&= \frac{{\bf{k}} \times {\bf{L}}}{\|\bf{L}\|} \\
&= \frac{{\bf{k}} \times ({\bf{x}}_B-{\bf{x}}_C)}{\|{\bf{x}}_B-{\bf{x}}_C\|}\\
&=\frac{-(y_B-y_C){\bf{i}}+(x_B-x_C){\bf{j}}}{\sqrt{(x_B-x_C)^2+(y_B-y_C)^2}}
\end{align}$$
$2.$ We can observe from the figure below that
$$d = ({\bf{x}}_B-{\bf{x}}_A)\cdot{\bf{n}}$$
so we have found the length of $AD$.
$3.$ The coordinates of $D$ will be
$$\begin{align}
{\bf{x}}_D &= {\bf{x}}_A + d{\bf{n}} \\
&= {\bf{x}}_A + [({\bf{x}}_B-{\bf{x}}_A)\cdot{\bf{n}}]{\bf{n}} \\
\\
&= {\bf{x}}_A + \left[({\bf{x}}_B-{\bf{x}}_A)\cdot\frac{{\bf{k}} \times ({\bf{x}}_B-{\bf{x}}_C)}{\|{\bf{x}}_B-{\bf{x}}_C\|}\right]\frac{{\bf{k}} \times ({\bf{x}}_B-{\bf{x}}_C)}{\|{\bf{x}}_B-{\bf{x}}_C\|} \\
\end{align}$$
and finally
$$\boxed{{\bf{x}}_D= {\bf{x}}_A + \frac{{\bf{k}} \cdot \left[ ({\bf{x}}_B-{\bf{x}}_C) \times ({\bf{x}}_B-{\bf{x}}_A) \right] \ }{\|{\bf{x}}_B-{\bf{x}}_C\|^2} {\bf{k}} \times ({\bf{x}}_B-{\bf{x}}_C)}$$
you can write the above vector equation as the following two scalar ones
$$\boxed{
\begin{align}
x_D &= x_A - \frac{\begin{vmatrix}
x_B-x_C & y_B-y_C \\
x_B-x_A & y_B-y_A \\
\end{vmatrix}}{(x_B-x_C)^2+(y_B-y_C)^2} (y_B-y_C) \\
\\
y_D &= y_A + \frac{\begin{vmatrix}
x_B-x_C & y_B-y_C \\
x_B-x_A & y_B-y_A
\end{vmatrix}}{(x_B-x_C)^2+(y_B-y_C)^2} (x_B-x_C)
\end{align}
}$$
Also, here is a MAPLE file that can help the reader of this post for the computation of final formulas mentioned here and the verification of the example mentioned in the question.
Best Answer
One simplification is to set \begin{align} a &= \text{xPos}_2 - \text{xPos}_1, \\ b &= \text{xVel}_2 - \text{xVel}_1, \\ c &= \text{yPos}_2 - \text{yPos}_1, \\ d &= \text{yVel}_2 - \text{yVel}_1. \end{align}
The vector of displacement from object $1$ to object $2$ at time $t$ is $(a + bt, c + dt)$, and the distance is $\sqrt{(a + bt)^2 + (c + dt)^2}$. Your function $f(t)$ becomes $$ f(t) = (a + bt)^2 + (c + dt)^2 = (b^2 + d^2)t^2 + 2(ab + cd)t + (a^2 + c^2). $$
In effect, what this does is to convert all your coordinates to the frame of reference of one of the two objects. The object whose frame of reference you choose is always at the point $(0,0)$ in that frame of reference, so all you have to track is the relative position of the second object, which travels along a straight line at constant speed.
Now it should be clear that $f(t)$ is the quadratic function $$ f(t) = At^2 + Bt + C \tag1$$ where $A = b^2 + d^2$, $B = 2(ab + cd)$, and $C = a^2 + c^2$. In the exceptional case where $b=d=0$, there is no relative motion, and $f(t)$ remains a constant $a^2 + c^2$. In any other case, you can find the exact minimum using derivatives, or simply recall the formula for the vertex of a parabola in this form, which tells us that the minimum occurs at $$ t_{\min} = -\frac{B}{2A} = -\frac{ab + cd}{b^2 + d^2}. $$ If you plug in known values of $\text{xPos}_1$, $\text{xPos}_2$, etc., this gives an exact solution for $t_{\min}$, which you can plug into equation $(1)$ to get the minimum distance without any guesswork or approximation.
Note that if you replaced $b_x$, $d_x$, $b_y$, and $d_y$ in your version of $f(t)$ with their equal multiples of $t$ (that is, do not replace expressions such as $t \cdot \text{xVel}_1$ with something "simpler"), expanded all the products, and combined terms with the same powers of $t$, you would eventually reduce $f(t)$ to something in the form $ f(t) = At^2 + Bt + C$ where $A$, $B$, and $C$ came out to the same values in terms of $\text{xPos}_1$, $\text{xPos}_2$, etc., as in the calculations above. It would just take more manipulation do to so.