This is a problem I can follow using summation notation as done in Merton (1972) but I'm having trouble finding the same solution using linear algebra methods.
The problem of finding the frontier of mean variance efficient portfolios has been summarized in matrix algebra form from Merton (1972) thus:
Let $\Omega$ be the matrix of covariances, and we assume it is positive definite and non-singullar.
\begin{align} (1) \qquad \min \omega '\Omega \omega \\ \text{s.t.} \quad & \omega ' \mu = \mu_p \\ & \omega ' \ \iota = 1 \\ \end{align}
where $\omega$ is a vector of porportions of each stock to be invested in the portfolio $\mu$ is the vector of expected returns for all assets in the portfolio and $\iota$ is the conforming vector of ones.
Using Lagrangian multipliers (1) can be written as
\begin{align} (2) \quad L = \omega ' \Omega \omega + \gamma_1(\mu_p - \omega ' \mu) + \gamma_2(1 - \omega ' \iota) \end{align}
where γ1 and γ2 are lagrange multipliers.
The standard FOCs for a critical point are:
\begin{align} (3) \qquad 2 \Omega \omega - \gamma_1 \mu - \gamma_2 \iota = 0 \end{align}
combining the restrictions above we find the solution
\begin{align} (4) \qquad \omega_p = g + h \mu_p \end{align}
Where
$$\begin{align} \quad g = \frac{1}{D}[B(\Omega^{-1} \iota) - A(\Omega^{-1} \mu)] \\ \text{and} \end{align}$$
$$\begin{align} \qquad & \quad h = \frac{1}{D}[C(\Omega^{-1} \mu) - A(\Omega^{-1} \iota)] \qquad \end{align}$$
and
$$\qquad A = \iota ' \Omega^{-1}\mu \qquad$$ $$\qquad B = \mu ' \Omega^{-1} \mu \qquad$$ $$\qquad C = \iota ' \Omega^{-1}\iota \qquad$$ $$\qquad D = BC - A^2 \qquad$$
Can anyone show me how we get from (3) to the solution (4)?
First solve (3) for the optimal weights by multiplying by the inverse $\frac{1}{2} \Omega^{-1}$ to obtain
$$\tag{*}\omega = \frac{\gamma_1}{2} \Omega^{-1}\mu + \frac{\gamma_2}{2} \Omega^{-1}\iota .$$
Enforcing the constaints provides two linear equations for the unknown multipliers $\gamma_1/2$ and $\gamma_2/2$,
$$\tag{**}\mu_p = \omega' \mu = (\Omega^{-1}\mu)'\mu \frac{\gamma_1}{2} + (\Omega^{-1}\iota)' \mu\frac{\gamma_2}{2}, \\ 1 = \omega' \iota = (\Omega^{-1}\mu)'\iota \frac{\gamma_1}{2} + (\Omega^{-1}\iota)' \iota\frac{\gamma_2}{2}.$$
Since $\Omega = \Omega'$ , we have $$(\Omega^{-1}\mu)' = \mu'(\Omega^{-1})' = \mu' (\Omega')^{-1} = \mu' \Omega^{-1}, \\ (\Omega^{-1}\iota)' = \iota'(\Omega^{-1})' = \iota' (\Omega')^{-1} = \iota' \Omega^{-1}. $$
Hence, (**) can be rewritten as
$$\mu'\Omega^{-1}\mu \frac{\gamma_1}{2} + \iota'\Omega^{-1}\mu\frac{\gamma_2}{2} = \mu_p, \\ \mu'\Omega^{-1}\iota \frac{\gamma_1}{2} + \iota'\Omega^{-1} \iota\frac{\gamma_2}{2}= 1.$$
With $A$, $B$ and $C$ as given, and noting that $A = \iota'\Omega^{-1}\mu =\mu'\Omega^{-1}\iota $ (since a scalar equals its transpose), the system of equations becomes
$$B \frac{\gamma_1}{2} + A \frac{\mu_2}{2} = \mu_p, \\ A\frac{\gamma_1}{2} + C \frac{\mu_2}{2} = 1.$$
From here, you should be able to solve for $\gamma_1/2$ and $\gamma_2/2$, substitute into (*) and obtain the optimal solution for $\omega$ which you write as $\omega_p$.