This is an interesting problem, and the anwers given may really help. But I miss a procedure to obtain definite answers.

Someone said, that is a trivial problem, middle-high-school level. I found that it is quite complicated if one tries to find general expressions. These get very quickly very complicated due to square-root terms and so on, so best is to treat it numerically.

Define two circles, a straight line and the slope of the tangent at circle 1

k1 = Sqrt[r1^2 - x^2];
k1n = -Sqrt[r1^2 - x^2];
k2 = Sqrt[r2^2 - (x - d)^2];
k2n = -Sqrt[r2^2 - (x - d)^2];
g = a x + b;
rul1 = a -> D[k1, x] /. x -> x1

The line must have the same value as the circle1 at position x1. This gives b

eq1 = k1 == g /. rul1 /. x -> x1
rul2 = Solve[eq1, b][[1, 1]]

and at x2 the line must have the same values as circle2

eq21 = k2 == g /. x -> x2 /. rul1 /. rul2
eq22 = k2n == g /. x -> x2 /. rul1 /. rul2

The slopes of the line at x1 and x2 must be equal to those of the circles (to get rid of the square-roots I square them)

eq3 = (D[k1, x] /. x -> x1)^2 == (D[k2, x] /. x -> x2)^2

Now define values for the problem

vals = {r1 -> 2.5, r2 -> .9, d -> 4.2};

and solve it

sol1 = NSolve[
  {eq21, eq3} /. vals,
  {x1, x2}
  ]
sol2 = NSolve[
  {eq22, eq3} /. vals,
  {x1, x2}
  ]    
sol1a = sol1 /. (u_ -> v_) :> v
sol2a = sol2 /. (u_ -> v_) :> v

Multiple solutions are possible, so define a test to find appropriate solutions

tst[y1_, y2_] := Module[{},
  aa1 = -(y1/Sqrt[r1^2 - y1^2]) /. vals;
  aa2 = -((y2 - d)/Sqrt[r2^2 - (y2 - d)^2]) /. vals;
  bb = r1^2/Sqrt[r1^2 - y1^2] /. vals;
  {aa1, aa2, bb}
  ]

and apply the test

pp1 = tst @@@ sol1a
param1 = Select[pp1, Abs[#[[1]] - #[[2]]] < .0001 &] // Flatten
pp2 = tst @@@ sol2a
param2 = Select[pp2, Abs[#[[1]] + #[[2]]] < .0001 &] // Flatten

and see what you have got

Plot[
 {
   k1, k1n,
   k2, k2n,
   param1[[1]] x + param1[[3]] /. vals,
   param2[[1]] x + param2[[3]] /. vals,
   -param1[[1]] x - param1[[3]] /. vals,
   -param2[[1]] x - param2[[3]] /. vals
   } /. vals,
 {x, 0, 6},
 AspectRatio -> Automatic
 ]

There is another somewhat shorter approach due to the fact that the straight line is given by the value of circle 1 and its tangent at x1

aa = -x1/Sqrt[r1^2 - x1^2]
bb = Sqrt[r1^2 - x1^2]
gg = aa (x - x1) + bb

At x2 the straight line must have the same value as circle 2

gla = (gg /. x -> x2) == Sqrt[r2^2 - (x2 - d)^2]

Solve for x2

lsgx2 = Flatten[ FullSimplify[Solve[gla, x2], And[r1 > 0, r2 > 0, d > 0]]]

x2 shall be unique, therefore the argument of the square-root has to vanish, giving possible values of x1

z1 = Cases[x2 /. lsgx2[[1]], Sqrt[x__] :> x, Infinity][[1]]
z2 = Cases[x2 /. lsgx2[[2]], Sqrt[x__] :> x, Infinity][[1]]
z1 - z2 // FullSimplify
lsgx1 = Solve[z1 == 0, x1]

with this x1 we can determine x2 (the solutions r1 and - r1 are obviously not appropriate)

lsgx2a = Flatten[FullSimplify[lsgx2 /. Take[lsgx1, -2]]] //  DeleteDuplicates

calculate the appropriate x1 and x2 and get the straight lines

lsgx1num = Take[lsgx1, -2] /. vals
lsgx2num = lsgx2a /. vals
sL1 = gg /. lsgx1num[[1]] /. vals
sL2 = gg /. lsgx1num[[2]] /. vals

and plot it

Plot[
 {
   k1, k1n,
   k2, k2n,
   sL1,
   sL2
   } /. vals,
 {x, 0, 6},
 AspectRatio -> Automatic
 ]

I tested only a few combinations of vals. The methods worked, but it is possible that there are combinations of circles where it may fail.

A few years ago I wrote this notebook which you may find helpful. It is in Spanish, but I have translated the relevant part for your question, if you go to button 3.

Here is a Manipulate[] to help you get started:

Manipulate[Graphics[{{Directive[AbsoluteThickness[2], ColorData[97, 1]],
                      Circle[c1, r], Circle[c2, r]}, 
                     With[{si = 2 r/EuclideanDistance[c1, c2]},
                          {Directive[AbsoluteThickness[3], ColorData[97, 4]], 
                           InfiniteLine[(c1 + c2)/2,
                                   {{Sqrt[1 - si^2], -si}, {si, Sqrt[1 - si^2]}} . (c2 - c1)], 
                           InfiniteLine[(c1 + c2)/2,
                                   {{Sqrt[1 - si^2], si}, {-si, Sqrt[1 - si^2]}} . (c2 - c1)]}],
                      {Directive[AbsolutePointSize[4], AbsolutePointSize[5]], Point[(c1 + c2)/2]}},
                     PlotRange -> 5],
            {{c1, {0, 0}}, Locator}, {{c2, {3, 2}}, Locator}, {{r, 1}, 1/20, 2}]

The important parts for you here would be the expressions within the InfiniteLine[] objects representing the internal tangents, since these can be used to give the equations of those lines, using e.g. RegionMember[InfiniteLine[(* stuff *)], {x, y}], which you can feed into Solve[] along with the equations for your circles to get the tangency points.

As to how those expressions were derived, these easily came out from the formulae in this Wiki article, which was greatly simplified by the fact that both of your circles have the same radii. (Another way to look at it is that the internal tangents both pass through the homothetic point of the two circles, which in the case of equal radii is the same as the midpoint of the two centers. From there, one can apply the Pythagorean theorem to reckon out the angle made by the internal tangents and the line joining the centers, recalling that one can draw a radius that is perpendicular to the tangent line of a circle.)

Good luck!