I have a few suggestions. You might try LibreCAD, which has a solver for common tangents to two ellipses, but I don't know anything about the API. The solver solves quartic equations, which is what you obtain if you naively try to find the common tangents for two ellipses.
If you want to roll your own: With a little bit of theory ("ranges of conics"), what you are asking for can be done with linear algebra (namely, finding inverses of 3x3 matrices) plus solving quadratics and one cubic equation. It goes like this:
You can express any conic (such as an ellipse) in the form of a matrix equation
[m00 m01 m02] [x]
[x,y,z] [m10 m11 m12] [y] = 0
[m20 m21 m22] [z]
where the matrix M is symmetric and [x,y,z] are homogeneous coordinates; just think z=1. We can write that equation in short form as X M X^T = 0 where X^T is the transpose of X.
Lines in the plane can be written in the form lx+my+nz=0 and so have "line coordinates" (l,m,n).
The set of tangent lines to the above conic can be expressed very simply in this notation. Let A be the inverse of the matrix M. The set of tangent lines to the conic is then
[a00 a01 a02] (l)
(l,m,n) [a10 a11 a12] (m) = 0
[a20 a21 a22] (n)
Now suppose we have a second conic with matrix N, and that N has inverse B. A common tangent will satisfy the above equation and the equation
[b00 b01 b02] (l)
(l,m,n) [b10 b11 b12] (m) = 0
[b20 b21 b22] (n)
In fact, we can scalar multiply the matrix in the latter equation by t and it will still hold:
[b00 b01 b02] (l)
(l,m,n) t [b10 b11 b12] (m) = 0
[b20 b21 b22] (n)
Adding the equation for the tangents to the first conic to the above equation for the second conic, we have the matrix equation L (A + tB) L^T = 0. So any common tangent to the two conics is a common tangent to every conic in the "range" A + tB.
Now for the big simplifying idea: we can find some very special conics in that range, "degenerate" conics, which are just pairs of points. Since the common tangents must pass through all conics, it follows that they must pass through the degenerate conics. But it is easy to find the lines that pass through degenerate conics, since such conics are just pairs of points!
You find the degenerate conics by solving the cubic equation det(A + tB) = 0 where det() is the determinant of 3x3 matrices. The cubic can be solved in closed form by Cardano's formula or a variation, or it can be solved numerically if that's what you need. Once you find the value(s) of t which make degenerate conics, you can factorize the equation L (A + tB) L^T = 0 into two linear factors. Each linear factor xl + ym + zn = 0 defines a point in homogeneous coordinates [x,y,z], or in Cartesian coordinates (x/z,y/z). You should get three point-pairs in this manner (six points). Taking lines through certain of those point-pairs will give you your four (at most) tangent lines.
Here is a simple example (where the centres of the two ellipses are both at the origin): find the common tangents to x^2+2y^2=3 and x^2+14y^2=7. The conic in matrix form are
[1 0 0] [x] [1 0 0] [x]
[x,y,z] [0 2 0] [y] = 0, [x,y,z] [0 14 0] [y] = 0
[0 0 -3] [z] [0 0 -7] [z]
The tangents are given by
[6 0 0] (l) [14 0 0] (l)
(l,m,n) [0 3 0] (m) = 0, (l,m,n) [ 0 1 0] (m) = 0
[0 0 -2] (n) [ 0 0 -2] (n)
Note I have multiplied the inverse matrices by a scalar just to make the entries integers rather than rational numbers. You don't have to do that but it makes the hand calculations easier. Multiplying the second matrix by an additional scalar t gives
[6+14t 0 0 ] (l)
(l,m,n) [0 3+t 0 ] (m) = 0
[0 0 -2-2t] (n)
The conic is degenerate when (6+14t)(3+t)(-2-2t)=0, i.e., when t=-3/7, -3, -1. When t=-3/7 we get
18/7 m^2 - 8/7 n^2 = 2/7 (9 m^2 - 4 n^2) = 2/7 (3m - 2n)(3m + 2n) = 0
This corresponds to the points with homogeneous coordinates [x,y,z] = [0,3,-2] and [0,3,2], in other words to points with Cartesian coordinates (0,-3/2) and (0,3/2).
When t=-3 we get -36l^2 + 4n^2 = (6l+2n)(-6l+2n) = 0, points [6,0,2] and [-6,0,2] or in Cartesian coordinates (3,0) and (-3,0). Finally, when t=1 we get -8l^2 + 2m^2 = 2(2l+m)(-2l+m) = 0 corresponding to points [2,1,0] and [-2,1,0] which are points at infinity.
Avoiding the points at infinity for now, just because they're a little more difficult to work with, we get four lines through the following pairs of points:
{(0,-3/2),(-3,0)}, {(0,-3/2),(3,0)}, {(0,3/2),(-3,0)}, {(0,3/2),(3,0)}
which give us the four common tangents to the two ellipses.
You can see from the picture that the common tangents also pass through the points at infinity [2,1,0] and [-2,1,0], i.e., that the pairs of parallel lines have slope 1/2 and -1/2.
Isn't that beautiful?
