@Axel Vogt 

It is not constant as

plots[odeplot](solA0,[theta,r(theta)],theta=varepsilon..5/1000);

shows. It is only a precision problem, abserr and relerr are used probably with some heuristic tests which are not enough when the solution is singular.

If you hope for some help, please state the complete problem, without asking to google for "research papers" with no references/links.

@tomleslie 

Is seems to be hard to obtain o prescribed precision.

In your worksheet one obtains the value of A:

A0 := 199.961318896644828500644760320

and the solution

Now, if we increase a bit the precision in dsolve without parameters for this A0:

solA0 := dsolve( {Ode1, r(varepsilon) = A0, (D(r))(varepsilon) = 0},
               numeric,
               range = varepsilon .. 5*(1/1000),
               abserr=1.0e-14,
               relerr=1.0e-13
             );

one obtains for solA0(5/1000);

i.e.  r(theta) has only 8 correct digits!

How many digits can we guarantee for A0?

@tomleslie

I do not think that the optimization is essential here; fsolve will do it after adjusting Digits.

The real problem is the accuracy of dsolve with parameters.

@Preben Alsholm 

I think we should have serious concerns about the precision because of the singularity.

For example fsolve (instead of NLPSolve) refuses to solve the equation even in the interval  199.96..199.97.

If we set Digits:=27 before fsolve ==>

199.961306129991433460330154

but the problem is how good is the approximation given by dsolve with parameters.

@acer 

I did not do it on purpose.

showstat(plottools:-`transform/object`);

simply does not include it.

Anyway I don't think that maplesoft will have something to lose, on the contrary! 

@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.