Why don't you use theta0(z)  ( or theta[0](z) ) instead of theta(z,0) ?

ode := diff(theta0(z), z, z)+theta0(z):
ibc:=theta0(0)=beta, D[1](theta0)(0)=0:
dsolve({ode,ibc});
                    theta0(z) = beta cos(z)

The result of fsolve could be NULL (if no solution found).

You should test this first, otherwise  NULL < 0   will produce the error.

pts:=[3,7],[6,5],[5,3],[1,2],[7,1]:
population:=5,2,2,1,3:

res:=seq(x-pts[i][1] <=  a||i,i=1..nops([pts])),
     seq(x-pts[i][1] >= -a||i,i=1..nops([pts])),
     seq(y-pts[i][2] <=  b||i,i=1..nops([pts])),
     seq(y-pts[i][2] >= -b||i,i=1..nops([pts])):
c:=add(population[i]*(a||i+b||i),i=1..nops([pts])):
s:=simplex[minimize](c,[res]):
eval(["x"=x,"y"=y,"min"=c],s);
                 ["x" = 5, "y" = 5, "min" = 51]

I always prefer a homemade version:

plot3d([u,(u^2-1)*cos(v)+1,(u^2-1)*sin(v)],u=1..2,v=0..2*Pi,
axes=normal,view=[-1..3,-4..6,-5..5],labels=[x,y,z],axes=boxed);

What is degrevlex?

Use e.g.

G:=Basis(K, 'tord', variables={r,u,v,w}); tord;

You could define the polynomial functions and use the composition operator @.

Or, better, use subs.
Notice first that the degrees of a and b must be 1.

a:=a1*x+a0;
                           a1 x + a0
b:=b1*x+b0;
                           b1 x + b0
c:=subs(x=a,b);
                      b1 (a1 x + a0) + b0
coeffs(collect(c-x,x),x):
sol:=solve([%],[a0,a1,b0,b1]):

# So, the general form of a,b is

{'a','b'} = eval({a,b},op(sol)); 

           {a,b} = {x/b1+a0, -a0*b1+b1*x}

# here a0,b1 are arbitrary constants, b1<>0.

Edit: actually c computes b o a, but the final result is of course the same. 

> min(extrema(x,{x^2 + x*y + y^2 = 2}));

It's a bug.

It seems that Apolonius works correctly only if all the 8 circles exist.