I've come across with the following proof of the Laplace expansion:
Let
$\Delta=\sum_{j=1}^n (-1)^{1+j} a_{1j}\bar M_j^1$
and
$\tilde{\Delta}= \sum_{j=1}^n (-1)^{i+j} a_{ij}\bar M_j^i$
We'll prove by induction that $\tilde{\Delta}=\Delta$.
For a $2 \times 2$ matrix, clearly, it holds:
$\tilde{\Delta_2}=\sum_{j=1}^n (-1)^{i+j} a_{ij}\bar M_j^i=-a_{21}a_{12}+a_{22}a_{11}=\Delta_2$
Suppose that for $n -1$ matrix $\tilde\Delta_{n-1}=\Delta_{n-1}$ is true. Then:
$$\tilde{\Delta_n}=\sum_{j=1}^n (-1)^{i+j} a_{ij}\bar M_j^i=\sum_{j=1}^n (-1)^{i+j} a_{ij}\left(\sum_{k<j}(-1)^{1+k}a_{1k}\bar M_{jk}^{i1}+\sum_{k>j}(-1)^k a_{1k}\bar M_{jk}^{i1}\right) $$
Gathering all the coefficients of $\bar M_{j_0k_0}^{i\,\,1}$
$$j_0>k_0\colon\; (-1)^{i+j_0}a_{ij_0}(-1)^{1+k_0}a_{1k_0}+(-1)^{i+k_0}a_{ik_0}(-1)^{j_0}a_{1j_0}=(-1)^{i+j_0+k_0+1}(a_{ij_0}a_{1k_0}-a_{ik_0}a_{1j_0})=(-1)^{i+j_0+k_0+1}M_{j_0k_0}^{1\,\,i}$$
$$j_0<k_0\colon\; (-1)^{i+j_0}a_{ij_0}(-1)^{k_0}a_{1k_0}+(-1)^{i+k_0}a_{ik_0}(-1)^{j_0+1}a_{1j_0}=(-1)^{i+j_0+k_0+1}(a_{ik_0}a_{1j_0}-a_{ij_0}a_{1k_0})=(-1)^{i+j_0+k_0+1}M_{j_0k_0}^{1\,\,i}$$
$$\tilde{\Delta_n}=\sum_{j\ne k}(-1)^{i+j+k+1}M_{jk}^{1i}\bar M_{jk}^{1i}$$
$$\Delta=\sum_{j=1}^n (-1)^{1+j} a_{1j}\bar M_j^1=\sum_{j=1}^n (-1)^{1+j} a_{1j}\left(\sum_{k<j}(-1)^{i+k-1}a_{ik}\bar M_{jk}^{1i}+\sum_{k>j}(-1)^{i+k-2}a_{ik}\bar M_{jk}^{1i}\right) $$
Gathering all the coefficients of $\bar M_{j_0k_0}^{i\,\,1}$
$$j_0>k_0\colon\; (-1)^{1+j_0}a_{1j_0}(-1)^{i+k_0-1}a_{ik_0}+(-1)^{1+k_0}a_{1k_0}(-1)^{i+j_0-2}a_{ij_0}=(-1)^{j_0+i+k_0}(a_{1j_0}a_{ik_0}-a_{1k_0}a_{ij_0})=(-1)^{i+j_0+k_0+1}M_{j_0k_0}^{1\,\,i}$$
$$j_0<k_0\colon\; (-1)^{1+j_0}a_{1j_0}(-1)^{i+k_0-2}a_{ik_0}+(-1)^{1+k_0}a_{1k_0}(-1)^{i+j_0-1}a_{ij_0}=(-1)^{k_0+i+j_0}(a_{1k_0}a_{ij_0}-a_{1j_0}a_{ik_0})=(-1)^{i+j_0+k_0+1}M_{j_0k_0}^{1\,\,i}$$
$${\Delta_n}=\sum_{j\ne k}(-1)^{i+j+k+1}M_{jk}^{1i}\bar M_{jk}^{1i}=\tilde{\Delta_n} ~~ \blacksquare $$
I just can't understand the part where they "gather coefficients":
$$j_0>k_0\colon\; (-1)^{i+j_0}a_{ij_0}(-1)^{1+k_0}a_{1k_0}+(-1)^{i+k_0}a_{ik_0}(-1)^{j_0}a_{1j_0}=\\ =(-1)^{i+j_0+k_0+1}(a_{ij_0}a_{1k_0}-a_{ik_0}a_{1j_0})=(-1)^{i+j_0+k_0+1}M_{j_0k_0}^{1\,\,i}$$
I'm okay with the $ (-1)^{i+j_0}a_{ij_0}(-1)^{1+k_0}a_{1k_0}$ part, which is simply what we get after plugging $j_0$ and $k_0$ into the previous formula of $\tilde{\Delta_n}$ (the sum of two sums) for a single and specific iteration. But where does the $$(-1)^{i+k_0}a_{ik_0}(-1)^{j_0}a_{1j_0}$$ come from? I just can't see anything like this before $\bar M_{j_0k_0}^{i\,\,1}$.
I'm trying to figure this out for many hours and I'm completely lost. I would appreciate any explanation!