I want to find $x_{11}$, $x_{12}$, $x_{21}$, $x_{22}$, $a$, and $b$ that solve the following system of equations \begin{align*} \begin{bmatrix}x_{11} & x_{12}\\x_{21} & x_{22}\end{bmatrix}\begin{bmatrix}1\\0\end{bmatrix}& = \begin{bmatrix}\gamma_{11} & \gamma_{12}\\\gamma_{21} & \gamma_{22}\end{bmatrix} \begin{bmatrix}x_{11} & x_{12}\\x_{21} & x_{22}\end{bmatrix}\begin{bmatrix}a\\0\end{bmatrix}, \\[2ex] \begin{bmatrix}x_{11} & x_{12}\\x_{21} & x_{22}\end{bmatrix}\begin{bmatrix}0\\1\end{bmatrix}& = \begin{bmatrix}\gamma_{11} & \gamma_{12}\\\gamma_{21} & \gamma_{22}\end{bmatrix} \begin{bmatrix}x_{11} & x_{12}\\x_{21} & x_{22}\end{bmatrix}\begin{bmatrix}0\\b\end{bmatrix}, \\[2ex] \begin{bmatrix}x_{11} & x_{12}\\x_{21} & x_{22}\end{bmatrix}\begin{bmatrix}\theta_{2}\\\theta_{1}\end{bmatrix}& = \begin{bmatrix}\alpha_{1}\\\alpha_{2}\end{bmatrix}. \end{align*} All parameters ($\gamma_{11}$, $\gamma_{12}$, $\gamma_{21}$, $\gamma_{22}$, $\theta_1$, $\theta_2$, $\alpha_1$, and $\alpha_2$) are non-zero and $\gamma_{11}\gamma_{22}-\gamma_{12}\gamma_{21}\neq0$.
Treating the system as 6 separete equations, I can express the solution as \begin{align*} x_{11} & =\frac{a\gamma_{12}}{1-a\gamma_{11}}x_{21},\\ x_{12} & =\frac{b\gamma_{12}}{1-b\gamma_{11}}x_{22},\\ x_{21} & =\frac{1}{\left( a-b\right) \theta_{2}}\left( \frac{\gamma _{11}\alpha_{2}-\gamma_{21}\alpha_{1}}{\gamma_{11}\gamma_{22}-\gamma _{12}\gamma_{21}}-b\alpha_{2}\right), \\ x_{22} & =\frac{1}{\left( b-a\right) \theta_{1}}\left( \frac{\gamma _{11}\alpha_{2}-\gamma_{21}\alpha_{1}}{\gamma_{11}\gamma_{22}-\gamma _{12}\gamma_{21}}-a\alpha_{2}\right), \end{align*} with $a$ and $b$ distinct roots of $$1-\left(\gamma_{11}+\gamma_{22}\right) z+\left(\gamma_{11}\gamma_{22}-\gamma_{12}\gamma_{21}\right) z^{2}=0.$$
I would like to represent the solution for $\begin{bmatrix}x_{11} & x_{12}\\x_{21} & x_{22}\end{bmatrix}$ in matrix form. Is this even possible? Any suggestion would be highly appreciated.
Let $G=[\gamma_{ij}]$ and $X=[x_{ij}].$
As the OP already figured out themselves, $a$ and $b$ are the roots of $1-\mathrm{tr}(G)\,z+\det(G)\,z^2=0.$ If we substitute $z=\lambda^{-1}$ and multiply with $\lambda^2,$ we get $\lambda^2-\mathrm{tr}(G)\,\lambda+\det(G)=0.$ This means that $a^{-1}$ and $b^{-1}$ are the eigenvalues of $G,$ and it might make sense to diagonalize $G.$ So let $$G = P \begin{pmatrix} a^{-1} & 0 \\ 0 & b^{-1}\end{pmatrix}P^{-1} = PDP^{-1} $$ with an invertible matrix $P,$ and let us see what we get. $$ X\begin{pmatrix} 1 \\ 0\end{pmatrix} =PDP^{-1}X\begin{pmatrix} a \\ 0\end{pmatrix} $$ or $$ P^{-1}X\begin{pmatrix} 1 \\ 0\end{pmatrix} =DP^{-1}Xa\begin{pmatrix} 1 \\ 0\end{pmatrix} =aDP^{-1}X\begin{pmatrix} 1 \\ 0\end{pmatrix} $$ which means $$ \left(I-aD\right)P^{-1}X\begin{pmatrix} 1 \\ 0\end{pmatrix} = 0 $$ This can easily be fulfilled iff the lower left element of $P^{-1}X$ is $0.$
From $$ X\begin{pmatrix} 0 \\ 1\end{pmatrix} =PDP^{-1}X\begin{pmatrix} 0 \\ b\end{pmatrix} $$ we get $$ \left(I-bD\right)P^{-1}X\begin{pmatrix} 0 \\ 1\end{pmatrix} = 0 $$ which means that the upper right element of $P^{-1}X$ must be $0.$ If we put all this together, we have $P^{-1}X=C$ with a diagonal matrix $C,$ which in turn means that $X=PC.$
Now we use the last of the three given equations to figure out the diagonal elements of $C.$ In other words, we figure out by which factor the initial choice of the eigenvectors of $G$ must be scaled to get the columns of $X.$
We want $X\theta = \alpha$ or $PC\theta = \alpha$ or $C\theta = P^{-1}\alpha$ with $\theta = (\theta_2,\theta_1)^T$ (sic!) and $\alpha=(\alpha_1,\alpha_2)^T.$
From this, we can conclude $$ C= \mathrm{diag}\left(P^{-1}\alpha\right) \, \mathrm{diag}\left(\theta_2^{-1},\,\theta_1^{-1}\right) $$ Therefore, we get $$ X= P \, \mathrm{diag}\left(P^{-1}\alpha\right) \, \mathrm{diag}\left(\theta_2^{-1},\,\theta_1^{-1}\right) $$