You could do as follows:

(beta - 1)*expand(remove(has,factor(e1),beta - 1));

Without going further into the details of your worksheet, there seems to be two problems with your metric M:
1.) Shold M[1,1] not have the opposite sign?
2.) In spherical coordinates, M[2,2] and M[3,3] cannot be put equal to 1. The Schwarzschild metric has M[2,2] = -r^2 and M[3,3] = -r^2*sin(theta)^2. The easiest way to access this metric is as follows:

Setup(metric = sc);

Do you not mean y = 0.6, rather than y = 0,6? With the former, there is no error raised.

Update: My objection above is erroneous as I had clearly misunderstood asympt (actually, I have never used it before). Looking at the help page, it is now clear to me that y = 0, 6 should be two arguments, the latter being the truncation order. So, please forget this answer of mine.

Without `or`, the statement served to if by seq(...) becomes

R[4] = R[1], R[4] = R[2], R[4] = R[3]

for i = 4, that is, an invalid boolean statement.

I guess the best way would be to make a function/procedure:

F := proc(H)
   ... the solving code ...
end proc:

that takes H as a parameter, i.e., to put the solving part of your program inside a function. After that you can either simply write

map(F,[100,125,140,160]);

or perform a loop like follows:

for H in [100,125,140,160] do F(H) end do;

PS: I assume that your question concerns the first part of your attached worksheet. As far as I can see, the variable SOL is not defined anywhere.

First of all what I think are some typos:
- Should the signs on the right-hand side of your expression for dR/dt not be the opposite ones?
- With the Hamiltonian you give, should dH/dq not equal 2*q + 2, rather than the opposite sign?

And now to the equation. It may be solved as follows:

de := diff(R(t,q,p),t) - (
   +diff(R(t,q,p),q)*p
   -diff(R(t,q,p),p)*(2*q + 2)
);
sol := pdsolve(de,R(t,q,p));

The solution is given in terms of an arbitrary function _F1, just as you find, although with different arguments. Having reached this point, the equation has been solved; from "And then I'm lost. How do I solve this PDE in maple? " you seem to think otherwise, or what?

PS: Without knowing why, I am presently uanble to copy-paste as an image the above solution; that is why it does not figure explicitly above.

This is probably only a workaround: Apply the following function to any of your expressions:

removeDaggers := (expr) -> eval(expr,map(x -> Dagger(x) = x,
   {d_[mu],M,mu5,mu,Pi[B],eta,i,aa[A] (* etc *)}   # Your hermitians
)):

I think there is something wrong with := in the vars assignment line. If I delete it, and reenter it standardly as : followed by =, the total derivative of Q(vars) comes out nonzero. I became aware of the problem by noting that on executing vars; on a line below your code, Maple returned just vars, i.e., vars unassigned.

There are some problems with your expression. First, it is ill-defined for B = 0, so the range B = 0..1 would not be expected to make sense. Second, the plot, if at all possible, would seem to need 4 dimensions to be rendered, 3 dimensions for the domain, and 1 dimension for the range/image. But perhaps what you want, I am just guessing now, is to plot several lines, for different values of A and B, in the same plot(?). If so, then you could do something like the following:

f := (A,B,x) -> A/B*x - A^2/B + B:
plots:-display([seq(seq(
   plot(f(A,B,x),x = 0..1)
,A = 1/4..1,1/4),B = 1/4..1,1/4)]);

Without being at all an expert in the theory of generalized functions, I know that products of distributions should be handled with care; in general, they are ill-defined. Let us consider your product, Dirac(x)*Heaviside(x), and abbreviate it to D(x)H(x), or just DH, for shorthand reasons. What is its action on a test function f = f(x)?

Approach I:

<DH,f> =    int(D(x)*H(x)*f(x)  ,x = -infinity..infinity)
       =(?) int(D(x)*[H(x)*f(x)],x = -infinity..infinity)
       =(?) H(0)*f(0)
       =    f(0)   (using your choice of the Heaviside)
       =    <D,f>,

so that DH = D, in fact.

Approach II (where D and H are interchanged):

<HD,f> =    int(H(x)*D(x)*f(x)  ,x = -infinity..infinity)
       =(?) int(H(x)*[D(x)*f(x)],x = -infinity..infinity)
       =(?) int(D(x)*f(x),x = 0..+infinity)
       =(?) lim int(D(x)*f(x),x = epsilon..infinity) for epsilon -> 0
       =    lim 0 for epsilon -> 0
       =    0,

so that here HD = 0, as Maple (perhaps) suggests. The question marks indicate, of course, equalities I am not completely certain of.

1.) Enter your equations:

eq1 := q[1,1] = -cos(theta[A])*sin(theta[B])*cos(theta[5]) - sin(theta[A])*sin(theta[5]);
eq2 := q[1,2] = +cos(theta[A])*sin(theta[B])*cos(theta[5]) + sin(theta[A])*cos(theta[5]);
# etc.

2.) Collect them, and apply the definitions of thetaA and thetaB:

eqs := eval({eq1,eq2},{
   theta[A] = theta[1] + theta[2],
   theta[B] = theta[3] + theta[4]
});

3.) Trying to solve (I am not all that confident that a solution for your set of coupled trigonometric equations can readily be found, but, of course, I may be mistaken):

sol := solve(eqs,{theta1,theta2,theta3,theta4,theta5});

This may perhaps not be the simplest or smartest solution, but you could introduce an auxiliary function as follows:

myD := (expr) -> `if`(expr::`^`,
   op(2,expr)*(op(1,expr))^(op(2,expr) - 1)*myD(op(1,expr)),
   D(expr)
):
h:= F-> expand(evalindets(myD(F)/F, specfunc(D), d-> op(d)*'h'(op(d)))):

Note that myD is now, of course, used in h. The definition of myD assumes, as you state, that the exponent is a constant; otherwise the formula does, of course, not hold.

The problem is the use of local P, which is not necessary. As far as I can see, the code may be written as follows:

Ps := (f,g) -> simplify(
   +diff(f,P[1])*diff(g,X)-diff(f,X)*diff(g,P[1])
   +diff(f,P[2])*diff(g,Y)-diff(f,Y)*diff(g,P[2])
):

As Cm is a Matrix, and not a function, write Cm[j,k] instead of Cm(j,k), etc.

Some less trivial examples, than having only zero entries, could be the following ones:

1.) Using the keyword 'symbol':

Matrix(4,3,symbol = a);

2.) Using the map (i,j) -> something:

Matrix(4,3,(i,j) -> a[i,j]),
Matrix(4,3,(i,j) -> a||i||j),
Matrix(4,3,(i,j) -> sin(i*j*x));