I seem to have stumbled upon a strange formula, and I'm not quite sure that I understand it. Consider the determinant standard Vandermonde matrix $${}^nD = \left|\begin{matrix} 1 & x_{0} & x^{2}_{0} & \cdots & x^{n}_{0} \\ 1 & x_{1} & x^{2}_{1} & \cdots & x^{n}_{1} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & x_{n} & x^{2}_{n} & \cdots & x^{n}_{n} \\ \end{matrix}\right|$$ In numerical analysis, we can often use this matrix (along with smaller, similar matrices) to generate the coefficients of the $n$-th degree interpolating polynomial for a general function $f$ passing through the set of points $\{(x_i,f(x_i))\}_{i = 0}^n$. Indeed, if we define the notation
$${}^nD_i(f) = \left|\begin{matrix} 1 & x_{0} & x^{2}_{0} & \cdots & x^{i-1}_0 & f(x_0) & x^{i + 1}_0 & \cdots & x^{n}_{0} \\ 1 & x_{1} & x^{2}_{1} & \cdots & x^{i-1}_1 & f(x_1) & x^{i + 1}_1 & \cdots & x^{n}_{1} \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots & \ddots &\vdots \\ 1 & x_{n} & x^{2}_{n} & \cdots & x^{i-1}_n & f(x_n) & x^{i + 1}_n & \cdots & x^{n}_{n} \\ \end{matrix}\right|$$
then the $n$-th order Lagrange approximation for the function $f$ can be obtained from the formula $$p(x) = f(x_0) + \frac{\det\left({}^1D_1(f)\right)}{\det(D)}(x - x_0) + \cdots + \frac{\det\left({}^nD_n(f)\right)}{\det(D)}(x - x_0)\cdots (x-x_{n-1}).$$
Now we can define an analogous notation wherein we replace not just the column corresponding to $x^i$ but also some other $x^k$ by
$${}^nD_{i,k}(f,g) = \left|\begin{matrix} 1 & x_{0} & x^{2}_{0} & \cdots & x^{i-1}_0 & f(x_0) & x^{i + 1}_0 & \cdots & x^{k-1}_0 & g(x_0) & x^{k + 1}_0 & \cdots & x^{n}_{0} \\ 1 & x_{1} & x^{2}_{1} & \cdots & x^{i-1}_1 & f(x_1) & x^{i + 1}_1 & \cdots & x^{k-1}_0 & g(x_1) & x^{k + 1}_0 & \cdots & x^{n}_{1} \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots & \ddots &\vdots & \vdots & \vdots & \ddots &\vdots\\ 1 & x_{n} & x^{2}_{n} & \cdots & x^{i-1}_n & f(x_n) & x^{i + 1}_n & \cdots & x^{k-1}_0 & g(x_n) & x^{k + 1}_0 & \cdots & x^{n}_{n} \\ \end{matrix}\right|$$
and it would seem that the following identity holds: $${}^nD_{i,k}(f,g)\:{}^nD = {}^nD_i(f)\:{}^nD_k(g) - {}^nD_i(g)\:{}^nD_k(f).$$
The first few cases are fairly easy to verify (some naive Gaussian elimination / symmetric polynomial tricks can handle up to $n = 5$ nicely, and I verified up to $n = 12$ with a computer), but I'm having trouble coming up with a general proof of this fact. If anyone has an idea as to how to show this, I would appreciate the insight. Thank you!
Hint
Say that the polynomial functions are given below:
$$ f(x) = \sum_{j=0}^{n} f_{j}x^{j} $$
Let $M$ be an identity matrix of size $n+1$ but with its $i$-th column replaced with $ \begin{bmatrix} f_{0} & f_{1} & \dots & f_{n} \end{bmatrix}^{T} $. Then we have the following relation:
$$ \begin{align} D_{i}(f)=DM \implies \det{\left(D_{i}(f)\right)} &= \det{\left(D\right)} \cdot \det{\left(M\right)} \\ &= \det{\left(D\right)} \cdot f_{i} \end{align} $$
--
Say that the second polynomial functions are given below:
$$ g(x) = \sum_{j=0}^{n} g_{j}x^{j} $$
Now let $M$ be an identity matrix of size $n+1$ but with its $i$-th and $k$-th columns replaced with $ \begin{bmatrix} f_{0} & f_{1} & \dots & f_{n} \end{bmatrix}^{T} $ and $ \begin{bmatrix} g_{0} & g_{1} & \dots & g_{n} \end{bmatrix}^{T} $ respectively. Then we have the following relation:
$$ \begin{align} D_{i,k}(f,g)=DM \implies \det{\left(D_{i,k}(f,g)\right)} &= \det{\left(D\right)} \cdot \det{\left(M\right)} \\ &= \det{\left(D\right)} \cdot \det{\left(\begin{bmatrix}f_{i} & f_{k} \\ g_{i} & g_{k}\end{bmatrix}\right)} \end{align} $$
--
Hope this can help you see why you have such identity. Next, what about $D_{i,j,k,...}(f,g,h,...)$?