In this specific case, you have, from your first equation:
$$ac=2$$
and, since $a,c\in\mathbb{Z}$, and $2>0$ you have the following possible pairs:
$$\begin{array}{|c|c|}
\hline
a & c\\
\hline
2 & 1\\
\hline
1 & 2\\
\hline
-1 & -2\\
\hline
-2 & -1\\
\hline
\end{array}$$
In the same way, since $b,d\in\mathbb{Z}$, we have that:
$$\begin{array}{|c|c|}
\hline
b & d\\
\hline
3 & -1\\
\hline
1 & -3\\
\hline
-1 & 3\\
\hline
-3 & 1\\
\hline
\end{array}$$
Now, due to the second equation:
$$ad+bc=-5$$
we can see that the solution is $(a,b,c,d)=(2,1,1,-3)$.
Now, in general, the problem of finding a factoriasation of $$p(x)=p_0+p_1x+p_2x^2,\ p_0,p_1,p_2\in\mathbb{Z}$$
in the form of:
$$E(Ax+B)(Cx+D)$$
with integers coefficients, boils down to finding rational roots of the equation $f(x)=0$, since, if $\frac{a}{b}$ and $\frac{c}{d}$ are the two roots of $f(x)=0$, then:
$$p(x)=p_2(x-\frac{a}{b})(x-\frac{c}{d})=\frac{p_2}{ad}(ax-b)(cx-d)$$
So, if $x=\frac{a}{b}$ is a rational root of $p(x)=0$ with $(a,b)=1$, we will have:
$$p\left(\frac{a}{b}\right)=0\Rightarrow p_0+p_1\frac{a}{b}+p_2\frac{a^2}{b^2}=0\Rightarrow b^2p_0=a(p_1b+p_2a)\Rightarrow b^2p_0|a\overset{(a,b)=1}{\Rightarrow}p_0|a$$
So, our first result is that $p_0$ should be a divisor of $a$, if $\frac{a}{b}$ is a root of $p(x)=0$. Now, we also note that - in case $a\neq0$, since this is a trivial case, due to the fact that $p_0=0$ and $p(x)=x(p_1+p_2x)$:
$$p_0+p_1\frac{a}{b}+p_2\frac{a^2}{b^2}=0\Leftrightarrow p_2+p_1\frac{b}{a}+p_0\frac{b^2}{a^2}=0$$
So, we can see that, as previously, $p_2|b$ is needed, as well.
In this way you can find several similar criteria.
Note: The conditions $p_0|a$ and $p_2|b$ are necessary conditions in order for $p$ to have rational roots, so it does not mean that $p$ has rational roots if $p_0|a$ and $p_2|b$, but that if $p$ has rational roots, then $p_0|a$ and $p_2|b$.
As those in the comments have said, the way to really solve such a problem is via Gaussian Elimination
In Python this is easily implemented either from first principles using lots of for loops and summations.
But an even easier (and faster!) solution is to use the NumPy package routine "numpy.linalg.solve".
You are quite correct that the matrix must be square and of full rank in order for you to be able to solve it.
[Aside note: A square matrix is full rank if and only if it's determinant is nonzero.]
Using numpy.linalg.solve(A,B) is neat as it will raise an error "LinAlgError" if A is singular or not square. The example given at the bottom of the documentation page is entirely self explanatory, so I will not insult your intelligence by going over it again!
So you could do something like this, but you'd need your chemist to already have written the thing up as a matrix or you could be kinder to them and write a little code to create the matrices from the equations they enter as plain text (e.g. https://stackoverflow.com/questions/45220032/how-to-balance-a-chemical-equation-in-python-2-7-using-matrices/45220319)
A = input("Enter the coefficents of your matrix: ")
B = input("Enter the ordinate values: ")
x = np.linalg.solve(A, B)
print x
Hope this helps, have fun and good luck! :-)
Best Answer
The next step is
$$\left(2^x\right)^2-10\cdot2^x+16=0$$ wich is a quadratic equation in $2^x$.
Then $2^x=5\pm3$ and $$x=1\text{ or }x=3.$$