@taizoon 

So all you know about theta(z,p)  is contained in    pde=0.
lambda(p) seems to be a given function.
The substitution p=0 is of no use.

You may consider theta = theta(z) as depending of the parameter p.

So,

ode := diff(theta(z),z,z) + lambda(p)^2*theta(z)+p*lambda(p)^2*(sin(theta(z))-theta(z));
s:=dsolve(ode=0,theta(z));


# it contains 2 solutions (given implicitely). Take e.g. the first one.

s[1];

That is all you can say about your theta.
keep in mind that the constants _C1 and _C2 depend on p.

@taizoon 

You did not post the whole problem and I cannot guess what you are trying to do.

@Markiyan Hirnyk 

The simple ideea is to use:

min{ |f(x)| + |g(x)| : x in X } =
min{ a + b : -a <= f(x) <= a,  -b <= g(x) <= b,  x in X, a in R, b in R}.

So, the standard simplex can be used if f,g are affine
(and provides an  _exact_ solution).

@Carl Love 

MultiInt:=proc(A,B)
local r:=NULL,x,Ax,Bx;
for x in {op(A)} intersect {op(B)} do
    Ax:=select(`in`,A,[x]);   # or     Ax:=select(u->evalb(u=x), A)
    Bx:=select(`in`,B,[x]);
    if nops(Ax)<nops(Bx) then r:=r,op(Ax) else r:=r,op(Bx) fi
od;
[r]
end;

Probably the efficiency is worse (it would be interesting to test) but you must admit that it is much easier to read/understand.

@Carl Love 

It does not work correctly
L([1,2,3,1,1,1,2],[1,1,2,2,2,2,2,2]) ;
should return [1,1,2,2]

Probably a "one line definition" for L is more subtle. Of course, a longger one is easy.

@Carl Love 

If multipicity is considered in the lists, one should use e.g.

L:=(A,B) -> select(u->member(u,B), A);

(actually this should be modified for the correct number of common multiplicity)

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