@Earl  I think that for your purposes a procedure will be useful that will allow you to easily generate a chain of tangent circles of any length. The parameters of the  Circles  procedure are: n - circle's number, R - radius of the original circle, r1 - radius of the larger of the small circles. The procedure returns the center of the n-th small circle:

restart;
Circles:=proc(n,R,r1)
option remember;
if n=1 then return [2*sqrt(R*r1),r1] else
simplify([Circles(n-1,R,r1)[1]*sqrt(R)/(sqrt(R)+sqrt(Circles(n-1,R,r1)[2])),Circles(n-1,R,r1)[2]*R/(sqrt(R)+sqrt(Circles(n-1,R,r1)[2]))^2]) fi;
end proc:

# Examples of use

S0:=seq(Circles(n,4,1), n=1..15);
S:=seq(plottools:-circle(Circles(n,4,1),Circles(n,4,1)[2]), n=1..15):
plots:-display(plottools:-circle([0,4],4),S, scaling=constrained, view=[0..5.5,0..2.5],size=[900,600]);

@Rouben Rostamian  Thank you for this great construction! I converted your comment into an answer and vote up.

@Earl  Using the general formulas from my answer above, it is easy to verify that regardless of the size of the smaller circles, the centers of these circles are always on the parabola  y=x^2/(4*R) , where  R  is the radius of the original large circle. This is clearly visible in the animation:

restart;
OneFrame:=proc(R,r1)
local O,dist,O1,O20,O2,O3,O4,O5,O6,k1,k2,k3,k4,P1,P2,P3,P4,P5,P6:
uses plots, plottools;
dist:=(A,B)->sqrt((A[1]-B[1])^2+(A[2]-B[2])^2):
O:=[0,R]: O1:=[x1,y1]: O2:=[x2,y2]: O3:=[x3,y3]: 
O1:=eval(O1,solve({dist(O,O1)=R+r1,y1=r1}, {x1,y1},explicit)[1]):
O20:=eval([x2,y2],simplify(solve({dist(O1,O2)=r1+r2, dist(O,O2)=R+r2,y2=r2},{x2,r2,y2}, explicit)[1])):
O2:=[2*sqrt(R)*sqrt(r1)*sqrt(R)/(sqrt(R)+sqrt(r1)),r1*(sqrt(R)/(sqrt(R)+sqrt(r1)))^2]:
k1:=simplify(sqrt(R)/(sqrt(R)+sqrt(O2[2])));
O3:=simplify([O2[1]*k1,O2[2]*k1^2]);
k2:=simplify(sqrt(R)/(sqrt(R)+sqrt(O3[2])));
O4:=[O3[1]*k2,O3[2]*k2^2]:
k3:=simplify(sqrt(R)/(sqrt(R)+sqrt(O4[2])));
O5:=[O4[1]*k3,O4[2]*k3^2]:
k4:=simplify(sqrt(R)/(sqrt(R)+sqrt(O5[2])));
O6:=[O5[1]*k4,O5[2]*k4^2]:
P1:=circle(O1,O1[2],thickness=2); P2:=circle(O2,O2[2],thickness=2); P3:=circle(O3,O3[2],thickness=2); P4:=circle(O4,O4[2],thickness=2); P5:=circle(O5,O5[2],thickness=2); P6:=circle(O6,O6[2],thickness=2);

display(circle(O,R,thickness=2,color=blue),P1,P2,P3,P4,P5,P6,plot(x^2/(4*R), x=0..10, color=green, thickness=2), pointplot([O1,O2,O3,O4,O5,O6],color=red,symbol=solidcircle,symbolsize=4),scaling=constrained, view=[-4..12,0..9] );
end proc:


plots:-animate(OneFrame,[4,r1],r1=0.4..4,frames=73);

@vv  Thank you for the solution!

@Rouben Rostamian  
Thank you for your interest and original solution!

Here is my short manual solution without trigonometry for an arbitrary right triangle:

We use the well-known formula for the radius of a circle inscribed in a right triangle  (a + b - c) / 2 . Further (see the figure) triangles  AEF = FKC , since they are similar, and the inscribed circles in them are the same. Therefore,  F  divides the segment  AC  in the ratio a : b . We get  AF = a*c / (a + b) . The similarity coefficient of triangle  AEF  with respect to triangle  ABC  is  k = c / (a+b) .
So we get the final answer  r = (a+b-c)*c / (2*(a+b))   (for the right triangle with sides 3, 4, 5 radius r = 5 / 7 )

Of course my and your answers for right triangles are the same:

is((a+b-c)*c/(a+b)/2=c/(a/(c-b)+2+b/(c-a))) assuming c^2=a^2+b^2;

                                       true
  

@minhthien2016  It is not difficult to do. I already wrote in one of the comments that each element of the list  L , i.e.  [x, y, z, u, v, w]  means a tetrahedron, in which  x, y, z  are the lengths of 3 edges coming out of one vertex,  u is opposite z ,  v is opposite y , w is opposite x  . 

@aaadoors  The condition you are writing about is necessary but not sufficient for the existence of such a tetrahedron. Read about the Cayley-Menger determinant in Wikipedia  https://en.wikipedia.org/wiki/Cayley%E2%80%93Menger_determinant  

Counterexample to your statement:
Bottom base -  (15,16,17).  Side faces -  (15,8,9), (17,9,10), (16,10,8)

@minhthien  You asked about the number of all solutions. The following code finds this number. This number is 4 times smaller than the length of the list  L , because each solution is represented in this list by 4 copies, depending on the choice of the tetrahedron vertex (the first 3 numbers in each sublist of  L  are the lengths of the 3 edges coming out of one vertex).

restart;
with(combinat):
primes := [seq(ithprime(i), i = 1 .. 20)];
S:=choose(primes, 6):
F:=(x,y,z,u,v,w)->`if`(LinearAlgebra:-Determinant(<0,1,1,1,1; 1,0,x^2,y^2,z^2; 1,x^2,0,u^2,v^2; 1,y^2,u^2,0,w^2; 1,z^2,v^2,w^2,0>)>0,true,false): k:=0:
for s in S do
for t in [seq(seq([s1[],s2[]], s2=permute(convert(s,set) minus convert(s1,set))), s1=choose(s, 3))] do
x,y,z,u,v,w:=t[];
if F(t[]) and min(x+y+z,x+v+u,w+y+u,w+v+z)>max(x+w,y+v,z+u) then k:=k+1; L[k]:=t;  fi;
od: od:
L:=convert(L, list): 
k/4;

@vv  Now I understand, thank you.

@vv  You still haven't explained why your code is skipping solutions (see my comment above "Missed solutions")

@vv  My solution  [3, 5, 13, 7, 11, 17]  is correct. Below is a detailed check. Tetrahedron ABCS with vertex S and base ABC. AS=3, BS=13, CS=5, AB=11, AC=7, BC=17. There is no triangle (5,13,7) here.

restart; # Check solution [3, 5, 13, 7, 11, 17]
f:=(X,Y)->sqrt((X[1]-Y[1])^2+(X[2]-Y[2])^2+(X[3]-Y[3])^2): # Distance beetween points X and Y
A:=[0,0,0]: B:=[x,y,0]: C:=[7,0,0]: 
sol1:=solve({f(A,B)=11,f(B,C)=17,y>0},explicit);
B:=eval([x,y,0],sol1);
S:=[x,y,z]:
sol2:=solve({f(A,S)=3,f(B,S)=13,f(C,S)=5,z>0},explicit);
S:=eval([x,y,z],sol2);
with(plottools): with(plots):
T:=textplot3d([[A[],"A"],[B[],"B"],[C[],"C"],[S[],"S"]], color=black, font=[times,16], align={left,above}):
display(curve([A,B,C,A,S,C]),line(S,B),T, color=red, thickness=2, scaling=constrained, axes=none);
f(A,B), f(A,C), f(B,C), f(A,S), f(B,S), f(C,S); # Check of lengths of sides

@vv   It seems your code is missing some solutions. In my answer for the first 10 primes, the first solution found is  [ 3, 5, 13, 7, 11, 17]   (22 solutions in total). But when I ran your code (slightly modified) to solve the same problem, I only got 5 solutions, and the solution above is missing:

restart;
primes := [seq(ithprime(i), i = 1 .. 10)]:
with(combinat):
comb:= choose(primes, 6):
IsTetra:=proc(L::list(posint))
  local d:=Matrix(4, {(1,2)=L[1], (1,3)=L[2], (2,3)=L[3], (1,4)=L[4], (2,4)=L[5], (3,4)=L[6]}, shape='symmetric');
  Matrix(3, (i,j)-> d[4,i]^2+d[4,j]^2-d[i,j]^2);  # Schoenberg
  LinearAlgebra:-IsDefinite(%, query = 'positive_definite')
end proc:
n:=0:
for s in comb do
  if IsTetra(s) then 
     n:=n+1; L[n]:=s;
  fi:
od:
L:=convert(L,list);
nops(L);

@nm  You have not responded to my answer in any way, probably you are not interested in this approach. I still think that the image of several stream lines together with several arrows on them very clearly characterizes the vector field for the corresponding system of differential equations. I have significantly edited my answer, adding several arrows.

@vv  Edges  x, y, z  come out from the same vertex, u is opposite z ,  v is opposite y , w is opposite x  . 

@Rouben Rostamian  Thanks for finding this example!  I tested the  Into_4_Equal_Areas  procedure on this triangle and it found both solutions. To get symbolic solutions I set in the code Digits=15 and then set in identify command to evalf[8] , because unfortunately fsolve doesn't work very well, when calculating with 15  digits it returns  0.166666631418626  instead  1/6=0.166666666666667 .

L:=[[[0,0],[1,0]],[[1,0],[1/2,8/9]],[[1/2,8/9],[0,0]]];
Sol:=Into_4_Equal_Areas(L);
identify(evalf[8](Sol));
plots:-display(<Pic(L,Sol[1]) | Pic(L,Sol[2])>);