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!
Best Answer
Start with a translation -- i.e., shift things to a slightly modified coordinate system so that $B$ is the origin $(0, 0)$:
$$ A \mapsto \left[A_x-B_x, A_y-B_y\right] = A'\\ B \mapsto \left[B_x-B_x, B_y-B_y\right] = B'\\ C \mapsto \left[C_x-B_x, C_y-B_y\right] = C'$$
Compute the angles formed by $A'$ and $C'$ with respect to the positive $x$ axis:
$$\theta_A = \tan^{-1} \frac{A_y-B_y}{A_x-B_x}$$ $$\theta_C = \tan^{-1} \frac{C_y-B_y}{C_x-B_x}$$
Compute the angle between $A'$ and $C'$:
$$\theta = \cos^{-1} \left( \frac{A'\cdot C'}{\|A'\|\|C'\|} \right)$$
Divide by two, and add this angle to the smaller of $\theta_A$ and $\theta_C$:
$$\phi = \min \left\{ \theta_A, \theta_C\right\}+\frac{\theta}{2}$$
Find a unit vector that has this angle in the shifted coordinate system, i.e.
$$\begin{align*} \tan \phi &= \frac{D_y}{D_x} \\ D_x^2+D_y^2 &= 1 \end{align*}$$
Alternatively, ignore the fact that it must be a unit vector, and just set $D_x = 1$ and compute $D_y$--it will still be on the same line, no matter what.
Finally, shift things back to original coordinates:
$$D = \left[ D_x + B_x, D_y + B_y \right].$$
Example:
$$A = [3,3] \\ B = [2,2] \\ C = [1,3]$$
Then,
$$A' = [3-2,3-2] = [1,1] \\ C' = [1-2,3-2] = [-1,1],$$ $$ \theta = \cos^{-1} \frac{ 1\cdot (-1)+1\cdot 1}{\sqrt{1^2+1^2}\sqrt{(-1)^2+1^2}} = \cos^{-1} 0 = \pi/2\ \;\;\; (\text{90 degrees}),$$ $$ \theta_A = \tan^{-1} 1 = \pi/4\ \;\;\; (\text{45 degrees})$$ $$ \theta_C = \tan^{-1} -1 = 3\pi/4\ \;\;\; (\text{135 degrees})$$
The smaller of $\theta_A$, $\theta_C$ is $\theta_A$, so finally,
$\phi = \frac{\pi}{4}+\frac{1}{2}\theta = \frac{\pi}{4}+\frac{\pi}{4} = \frac{\pi}{2}$
as expected.
If $\phi = \pi/2$, then this is a vector pointing along the y-axis, so let $D' = (0,1)$, which leads to $D = D'+B = (0+2,1+2) = (2,3).$
Plot $A$, $B$, $C$ and D$ and you will see that this is the desired result.