This is strictly speaking not answering your question, but you may do as follows:

number := (x::posint) -> floor(log10(x) + 1):

You can use delayed evaluation:

'A' = A,'B' = B,'C' = C,'E' = E;

There is no need to explicitly define either the spacetime dimension or the metric, as the ones you choose are the default ones. I think you code can thus be simplified to

with(Physics):
Define(nn[mu] = Vector([1,0,0,0]),k[mu],quiet):
SumOverRepeatedIndices(nn[mu]*k[~mu]);

Somewhat a variation on Kitonum's theme, using select and remove:

s,r := selectremove(z -> has(z,x) and has(z,y),expr):
``(simplify(select(has,r*expand(s),x)))*
``(simplify(remove(has,r*expand(s),x)));

As Joe Riel has already pointed out, the result is correct. What you need to do is simply to interchange the two print statements, for the statement corresponding to the then-clause is by definition the one being executed if the boolean expression (in your case a > b) evaluates to true.

By defining an auxiliary type, T say, you can transform the expression using evalindets. Below an example with more than one parameter which is not to be included in the combine operation:

expr :=
   -sin(alpha)*(sin(theta1)*cos(theta)-cos(theta1)*sin(theta))
   -cos(beta) *(sin(theta1)*sin(theta)+cos(theta1)*cos(theta));
`type/T` := x -> not has(x,{alpha,beta}):   # alpha and beta not to be included
evalindets(expr,T,combine);

Perhaps the following is of interest to you:

sol := solve({x+y+z=1,x-y+3*z=7}):
lhs~(select(lhs = rhs,sol));   # The set of free parameters in sol

Update/warning: The output below results when the dimension and signature, respectively, are chosen as 3 and `+++`. This, of course, makes little sense, as the dimensionality is four; yesterday I obviously forgot to delete or uncomment these lines in my worksheet after having played around a bit. The difference in output disappears (including the issue of covariant/contravariant pairs) if these assignments are not made, and thus this post of mine seems to become quite vacuous. I became aware of this a moment ago when updating the Physics-package due to Edgardo's post, in order to check the triply contracted expression of illuminates.

I would prefer to have contracted pairs of indices to be properly covariant/contravariant pairs. And, in fact, it seems to make a difference [your expression is expr, mine with the properly matching indices is expr_]:

expr := LeviCivita[4,sigma,lambda,rho]*LeviCivita[4,xi,eta,mu]*g_[ rho, mu]*qp[ sigma]*q[ lambda]*qp[ xi]*q[ eta];
expr_:= LeviCivita[4,sigma,lambda,rho]*LeviCivita[4,xi,eta,mu]*g_[~rho,~mu]*qp[~sigma]*q[~lambda]*qp[~xi]*q[~eta];
Simplify(expr );
Simplify(expr_);
SumOverRepeatedIndices(expr );
SumOverRepeatedIndices(expr_);

But this state of affairs is not reassuring, I think. Checking the two expressions with Physics:-Check does not raise any warning:

Check(expr );
Check(expr_);

`The repeated and free indices in the given expression check ok.`
`The repeated and free indices in the given expression check ok.`

I think the following will do the trick:

with(Typesetting):
Settings(typesetdot = true):

But, coming to think of it, perhaps it is not available in Maple 18?

Here is one way of doing it:

de  := diff(1/r(phi),phi$2) + 1/r(phi) - G*M/h^2;
sol := dsolve({de,r(0)= 2/3,D(r)(0) = 0});
assign(sol):
plots:-polarplot([eval(r(phi),{
   M = 1,
   G = 1,
   h = 1
}),phi,phi = 0..2*Pi]);

If I understand your flowchart correctly, I would do something as follows:

for i from 1 by 1 while i <= 330 and max(Q) > 0 do
   do_function
end do;
if i <= 330 then print(Solution)   # While loop terminated due to max(Q) <= 0
elif max(Q) > 0 then print(NFS)
else print(Cycling)
end if:

PS: Note that the above code will only run if proper values have been assigned to Q, etc.

This is probably not what you have in mind, but I would like to mention the socalled jet notation:

PDEtools:-ToJet(di,y(t));

Using applyop, as you suggest, in combination with numer and denom:

applyop(x -> expand(x/denom(a)),3,numer(a));

Applying some assumptions [I have no idea whether these assumptions are actually valid in your case]:

simplify(original) assuming beta::positive,L__1::positive;

If P as well is assumed positive, then the following expression results:

A somewhat simplistic approach:

P := product(exp(1/i),i = 1..n);
for j from 1 to 60 do
   expr := evalf(eval(P,n = j));
   if expr > 100 then print(j,expr) end if;
end do:

It shows that the smallest value is n = 56.