Provided only with the radii and centers of two circles. Set up a system of equations for the points of tangency such that it can be solved by an iterative, numerical, method.
I have found this page http://mathworld.wolfram.com/Circle-CircleTangents.html. However, I am not sure how to translate those four equations to one system of equations (one matrix).
$$ (t_2-x_2)\cdot(t_2-t_1)=0\\ (t_1-x_1)\cdot(t_2-t_1)=0\\ \lvert t_1 -x_1 \rvert^2 = r_1^2\\ \lvert t_2 -x_2 \rvert^2 = r_2^2 $$
This might not be what you asked for, since I'm not using numerical methods, only simple square roots. I'm using a completely different approach than you'd likely have considered, though. I consider it rather elegant, so I'm posting this since it might be what you (or someone else reading your question) wants, if you only knew this was possible.
I'm going to address this question using techniques from Lie geometry, using the problem of Appollonius as the central building block.
You can represent a circle as a vector $(x, y, \pm r, x^2+y^2-r^2, 1)^T$ or any multiple thereof (homogeneous coordinates). Any circle $c$ will satisfy $c^TLc=0$, with the matrix $L$ (describing the Lie quadric) given by
$$L=\begin{pmatrix} -2 & 0 & 0 & 0 & 0 \\ 0 & -2 & 0 & 0 & 0 \\ 0 & 0 & 2 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 \\ 0 & 0 & 0 & 1 & 0 \end{pmatrix}$$
(I'm using a different basis than Wikipedia and most of the literature does, but I find it convenient.) Points are special circles of radius $0$, i.e. vectors of the form $(x,y,0,x^2+y^2,1)$ or multiples thereof. Lines are special circles of radius infinity and with a center at infinity as well. In homogeneous coordinates, this results in the last coordinate becoming zero. More precisely, a line $ax+by+c=0$ is encoded using the vector $(a,b,\pm\sqrt{a^2+b^2},-2c,0)^T$.
Two circles $a,b$ are in oriented contact to one another if they satisfy $a^TLb=0$. Oriented contact means that they touch one another and their orientation agrees. You can imagine orientation as an arrow drawn anywhere along the circle. If at the touching points the arrows agree in direction, you have oriented contact. So all circles are oriented, and the lines as well. That's the reason I had $\pm$ for the radius up there: change of sign switches orientation. Points have no radius and therefore no orientation.
The problem of appollonius is this: given three circles $a,b,c$, find a circle $d$ which touches all of these. So you want $a^TLd=b^TLd=c^TLd=d^TLd=0$. The first three of these can be encoded in a linear system of equations. This will have a two-dimensional space of solution. The last is a quadratic equation which will in general select two one-dimensional spaces from this set of solutions. Picking any non-zero representative from each of them gives you homogeneous coordinates for the touching circles. So even after you fix all the orientations (except perhaps for the result), you still get two possible answers. Changing signs of the input will change the output, but changing all signs in the input simply changes all signs in the output but leaves everything else in place.
In your case, you only have two circles, not three, but you also have the extra requirement that the touching “circles” should in fact be lines, with a zero in the last coordinate. So you can write $L$ times the two circles as two equations, and a third equation (representing the “point at infinity”) for the last equation. If you want a marix representation, I'd write it like this:
$$\begin{pmatrix} -2x_1 & -2y_1 & +2r_1 & 1 & x_1^2+y_1^2-r_1^2 \\ -2x_2 & -2y_2 & \pm2r_2 & 1 & x_2^2+y_2^2-r_2^2 \\ 0 & 0 & 0 & 0 & 1 \end{pmatrix}\cdot v=0$$
This is an underdetermined system of linear equations, which will in general have a two-dimensional solution space. Or in other words, the matrix has a two-dimensional right kernel. Name that kernel $\operatorname{span}(v_1,v_2)=\{\lambda v_1+\mu v_2\mid \lambda,\mu\in\mathbb R\}$. Find a linear combination which satisfies
$$(\lambda v_1+\mu v_2)^T\cdot L\cdot(\lambda v_1+\mu v_2) = \lambda^2(v_1^TLv_1)+2\lambda\mu(v_1^TLv_2)+\mu^2(v_2^TLv_2) = 0$$
with the matrix $L$ as given above. This is still a homogeneous equation. In general you can choose one of the two coefficients arbitrarily, for the other you will have two solutions. Each of the resulting vectors has the form
$$(\lambda v_1+\mu v_2)=(a,b,\pm\sqrt{a^2+b^2},-2c,0)^T$$
and corresponds to a line with the equation $ax+by+c=0$ which is tangent to your circles. As you saw my original linear system of equations had one sign given as $\pm$. The positive sign will lead to the outer tangents (both circles oriented the same way), the negative one to the inner tangents (circles oriented in opposite ways).
So now you have the lines, but you want the points of tangency. If you have two circles $a,b$ in oriented contact (i.e. $a^TLb=0$), then the contact element, i.e. the set of all circles touching these two in the same point with the same orientation, is simply the span of the two vectors. So take one circle and one tangent, and from their span choose the one element which is a point, i.e. has radius zero.
More precisely, if you have the circle $(x,y,\pm r,x^2+y^2-r^2,1)^T$ and the line $(a,b,\pm\sqrt{a^2+b^2},-2c,0)^T$, multiply the former by the third coordinate of the latter, and the latter by the third coordinate of the former, then subtract these two. You get something point-like, which you can divide by its last component to obtain a representative whose first two coordinates are the coordinates of the point of tangency.
The above contains the general approach I'd choose to any problem which can be addressed using Lie geometry. It's essentially a recipe for solving Apollonius' problem or anything closely related to this. In your particular case, the last line of the system of linear equations forces the last coordinate of any solution $v$ to be zero. This in turn means that the last column of the first two rows is irrelevant. When computing $\lambda$ and $\mu$ you could likewise ignore the third and fourth coordinate, and concentrate on the first three. But I'd like to show you where this came from, and how general the tool is.
While my other answer uses a completely different approach, this one here is closer to your description of the problem. Start with the equations you gave:
\begin{align*} 0 &= (t_2-x_2)\cdot(t_2-t_1) = t_2\cdot t_2 - x_2\cdot t_2 - t_2\cdot t_1+x_2\cdot t_1 \\ 0 &= (t_1-x_1)\cdot(t_2-t_1) = t_1\cdot t_2 - x_1\cdot t_2 -t_1\cdot t_1 +x_1\cdot t_1 \\ r_1^2 &= (t_1-x_1)\cdot(t_1-x_1) = t_1\cdot t_1 - 2\,t_1\cdot x_1 + x_1\cdot x_1 \\ r_2^2 &= (t_2-x_2)\cdot(t_2-x_2) = t_2\cdot t_2 - 2\,t_2\cdot x_2 + x_2\cdot x_2 \end{align*}
You asked for a matrix, which is typically used for linear equations. So you could write this as something like this:
$$ \begin{pmatrix} x_{2x} & x_{2y} & -x_{2x} & -x_{2y} & 0 & 1 & -1 \\ x_{1x} & x_{1y} & -x_{1x} & -x_{1y} & -1 & 0 & 1 \\ -2x_{1x} & -2x_{1y} & 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & -2x_{2x} & -2x_{2y} & 0 & 1 & 0 \end{pmatrix} \cdot\begin{pmatrix} t_{1x} \\ t_{1y} \\ t_{2x} \\ t_{2y} \\ \lVert t_1\rVert^2 \\ \lVert t_2\rVert^2 \\ t_1\cdot t_2 \end{pmatrix} =\begin{pmatrix} 0 \\ 0 \\ r_1^2 - \lVert x_1\rVert^2 \\ r_2^2 - \lVert x_2\rVert^2 \end{pmatrix} $$
If you ignore the inter-dependence between the entries in that seven-element vector, you have four equations defining seven variables, which in general leads to a three-dimensional space of solutions. So you could talk about linear combinations of these three vectors, and search for the ones satisfying the extra three components. But this doesn't simplify the problem that much: instead of looking for four variables determined by four quadratic equations, you'd now be searching for three variables determined by three quadratic equations. Therefore, you might as well solve your original problem directly, without ever writing it down as a a matrix.
So how would you feed all of this to numerical methods? I'd start by turning the four equations into a single one. Take the difference between left hand side and right hand side, and square that. Sum up the four squared differences. If you have a solution, that sum will be zero. As long as one equation is not satisfied and all your parameters are real, the sum will be positive. So now optimize over a four-dimensional parameter space ($x$ and $y$ coordinates of $t_1$ and $t_2$) with the objective of minimizing the sum of four squares. No guarantees that this will be well-behaved, but it's a starting point. Personally I'd stick to my other approach, though.
