Thank you for the info. I know that Geogebra is also popular. They include CASes, but it would be nice take advantage of a full featured system such as Maple.

It is not dificult to program it.

Here is a not optimal one where the circles are represented as lists [x,y,r].

myapo:=proc(c1,c2,c3)
local i1,i2,i3,x,y,r,s,u,c,apo:=NULL;
for i1 in [-1,1] do for i2 in [-1,1] do for i3 in [-1,1] do
s:=solve({ (x-c1[1])^2+(y-c1[2])^2 = (r*i1 + c1[3])^2,
       (x-c2[1])^2+(y-c2[2])^2 = (r*i2 + c2[3])^2,
       (x-c3[1])^2+(y-c3[2])^2 = (r*i3 + c3[3])^2 },
       {x,y,r});
for u in [s] do
    c:=evalf([eval(x,u),eval(y,u),eval(r,u)]);
    if not(type(c,list(realcons))) or (c[3]<0) then next fi;
    apo:=apo,c od;
od od od;
print(cat("number of circles:",nops([apo])));
apo
end:

To test it with your data:

c1,c2,c3:=[0,0,5],[5,4,2],[13,0,3]:

c8:=myapo(c1,c2,c3);
ccc:=[seq(circle([u[1],u[2]],u[3],color=red),u=[c1,c2,c3]), seq(circle([u[1],u[2]],u[3],color=black),u=[c8])]:
with(plottools):with(plots):
display(ccc);

@rnoxy 

You may use

try
A := Apollonius( c1, c2, c3 )
catch:   A:= MyApolonius(c1,c2,c3)
end try;

which catches any exception.

In this caze z=99/100 is not so close to 1 when the precision is 50 to explain the bug. Simply adding 10000 terms of the series provides 40 digits.

As a workaround:

> convert(hypergeom([2/3, 2/3, 2/3], [-1/3, 4/3], 99/100),MeijerG):
> evalf[50](%);

   -138.99623131713330917549791383444227380249824260751

   +8.3550383462461792073697301492096375870611508951167*10^(-47)*I