Example:

r:=rand(1..6):
L:=['r'() $ 30];
ListTools:-Occurrences(5, L);

Use the functional assignment:

f₀ := eta -> 1 - exp(-eta);

1. A function  f  can be applied only to one object or to a list (or to a set or to array) of objects by elementwise operation.

2. A argument (arguments) of a function  f  should be inside of parenthesis.

f([seq(i, i = 1 .. 2)][]);                      
f(add(i, i = 1 .. 2));

Three simple procedures  IsReflexive, IsSymmetric, IsTransitive  solves your problems.

IsReflexive:=proc(X::set, Y::set(list))
local k, y;
k:=0;
for y in Y do
if y[1]=y[2] then k:=k+1 fi;
od;
if nops(X)=k then true else false fi;
end proc:

Examples of use:

IsReflexive({a,b,c}, {[a,a],[b,b],[c,c],[a,b],[b,a]});
IsReflexive({a,b,c}, {[a,a],[b,b],[a,b],[b,a]});
                              true
                             false

IsSymmetric:=proc(X::set, Y::set(list))
local k, y;
k:=0;
for y in Y do
if [y[2],y[1]] in Y then k:=k+1 fi;
od;
if k=nops(Y) then true else false fi;
end proc:

Examples of use:
IsSymmetric({a,b,c},{[a,a],[b,b],[a,b],[b,a]});
IsSymmetric({a,b,c},{[a,a],[b,b],[a,b],[b,a],[a,c]});
                              true
                             false

IsTransitive:=proc(X::set, Y::set(list))
local y1, y2;
for y1 in Y do
for y2 in Y do
if y1[2]=y2[1] and not ([y1[1],y2[2]] in Y) then return false fi;
od; od;
true;
end proc:

Examples of use:
IsTransitive({a,b,c}, {[a,a],[b,b],[a,b],[b,c]});
IsTransitive({a,b,c}, {[a,a],[b,b],[a,b],[a,c]});
IsTransitive({a,b,c}, {[a,a],[b,b],[a,b],[b,c],[a,c]});
                             false
                              true
                              true
 

In Maple 2015 - 2017 there are no any changes:

restart;
(9*b^2 - r^2 - a^2*cos(theta)^2)^2;
                                           

You probably work with an older version of Maple. Try this option:

restart;
``(9*b^2 - r^2 - a^2*cos(theta)^2)^2;

@Daniel Skoog  A nice application (thumb up)!

Here's an alternative solution. For greater clarity, I made the cuboid opaque and increased the number of frames by 20 times. To create the animation,  more convenient  plots:-animate  command is used. The same plot shows the height change  h  and the value of the corresponding volume  V  (the height changes in 0.01 increments). Using the animation panel for step-by-step displaying of frames, we can see that the largest volume will be approximately equal to 1.19 for a height equal to 0.67  (I took the height of the pyramid to 2):

restart;
with(plots): with(plottools): 
Digits:=6:
Y:=[[-1,-1,0],[1,-1,0],[1,1,0],[-1,1,0],[0,0,2]]:

Pyramid:=display(curve([seq(y,  y=Y[1..-2]), Y[1]], thickness=3, color="Blue"), seq(line(y,Y[-1], color="Blue", thickness=3), y=Y[1..-2])):
 
animate(display, [cuboid([(2-a)/2,(2-a)/2,0],[-(2-a)/2,-(2-a)/2,a], color="Yellow"), textplot3d([[-0.5,1.4,2, V=a*(2-a)^2], [-0.5,1.4,2.2, h=a]], font=[times,18])], background=Pyramid, a=0..2, frames=201, axes=none, scaling=constrained, orientation=[60,70], paraminfo=false, lightmodel=light1);

Here is a screenshot of the animation frame, which shows the cuboid of the largest volume:

maximize(h*(2-h)^2, h=0..2, location);  # Exact calculation of the largest cuboid volume
evalf(%);

                              1.18519, {[{h = .666667}, 1.18519]}
 

I do not know the reason why  applyrule  does not work, because the documentation for this command is very poor. Use  evalindets  command instead:

restart;
eq1 := a[n+2]+4*a[n+1]+10*a[n]-5*n-1;
evalindets(eq1, `indexed`, t->op(t)*p+q+b[op(t)]);

Try

D[1](u)(0,t)=0, D[1](u)(L,t)=0

Addition.  I think that your teacher's proposal is incorrect in any version of Maple. Here is an example:

restart;
u:=(x,t)->x^2+3*x*t+t^4:
diff(u(0,t),x), diff(u(L,t),x);  # An incorrect way
# Two valid options:
D[1](u)(0,t), D[1](u)(L,t);   # The first way (the best one)
eval(diff(u(x,t),x), x=0), eval(diff(u(x,t),x), x=L);   # The second way

The notation that you are proposing is also incorrect in Maple, because derivative must first be calculated at an arbitrary point. Here is an valid way:

restart;
Du:=diff(u(x,t),x);
eval(Du, x=0);

Formally to extract the sequence of used arithmetic operators from  an expression you can as follows:

a:=3*x^3-5*x^2+3*y:
S:=convert(a,string):
Operators:={"+", "-", "*", "^"}: 
StringTools:-Select(s->is(s in Operators), S);
seq(s, s=%);
ListTools:-Collect([%]);

Edit.

Give these points names  A, B, C, E  that mean their coordinates:

curve := y^2 = x^3 - 43*x + 166;
Points:=[[3,8],[-5,16],[11,32],[3,-8]]:
A, B, C, E:= op(Points);
plots:-display([
plot(+sqrt(rhs(curve)),x = -10..12),
plot(-sqrt(rhs(curve)),x = -10..12),
plots:-pointplot([A, B, C, E],symbol = solidbox)
]);

To label these points on the plot, use plots:-textplot command.

Addition. After you give names to these points, when using plots:-textplot command, you do not need to retype their coordinates (you can just refer to their names):

restart;
curve := y^2 = x^3 - 43*x + 166:
Points:=[[3,8],[-5,16],[11,32],[3,-8]]:
Labels:=["A","B","C","E"]:
plots:-display([
   plot(+sqrt(rhs(curve)),x = -10..12),
   plot(-sqrt(rhs(curve)),x = -10..12),
   plots:-pointplot(Points, symbol = solidbox, color=[green,red,    blue,yellow], symbolsize=17), 
   plots:-textplot([seq([op(Points[i]),Labels[i]], i=1..4)],  align=above, font=[roman,16])
]);
 

Edit.

If you want Maple not to open the parentheses, then do this:

z:=-x-y*I;
(-1)*``(-z);
                            z := -x - y I
                             - (x + y I)

or immediately type

z:=-``(x+y*I);
                           z := - (x + y I)

If later you want Maple to open parentheses, then use  expand command:

z:=-``(x+y*I);
expand(z);
                                 -x - y I

This way you can always сarry the desired factor out of brackets of some  expression, for example:

P:=-3/2*x^2-x+5/2:
-1/2*``(P/(-1/2));
                                        -1/2*(3*x^2+2*x-5)

Edit.

The domain of definition of your function is a semicircle. Use the direct plotting in this domain:

plot3d(4*x*y^2-x^2, x=0..1, y=-sqrt(1-x^2)..sqrt(1-x^2), scaling=constrained, numpoints=10000, axes=normal);

Another option is the using of the polar coordinates:

plot3d(eval([x,y,4*x*y^2-x^2],[x=r*cos(t),y=r*sin(t)]), r=0..1, t=-Pi/2..Pi/2, scaling=constrained, axes=normal);

VerticesOfCuboid:=proc(V1::list, V2::list)
[V1, seq(subsop(i=V2[i], V1), i=1..3), V2, seq(subsop(j=V1[j], V2), j=1..3)];
end proc: 

Example of use:

VerticesOfCuboid([1,1,2], [3,4,5]);

               [[1, 1, 2], [3, 1, 2], [1, 4, 2], [1, 1, 5], [3, 4, 5], [1, 4, 5], [3, 1, 5], [3, 4, 2]]

Addition. In geom3d you can specify a cuboid by  parallelepiped  command and find its parameters by other commands (detail, vertices, faces  and so on).

There are several ways to do this. Here are two ones:

restart;
f := unapply(x*(8*x^2+5*x+cos(y)), x, y);
p := solve({diff(f(x, y), x) = 0, diff(f(x, y), y) = 0}, {x, y});
eval([x,y], p[1]);  # The first way
assign(p[1]);  # The second way 
x, y; 

I understood  this  to mean that  P  is an operator in the functional space and therefore we can do so:

P:=u->(i,j)->u(i+1,j);

Example of use:

u:=(i,j)->i^2+j^2;
v:=P(u);  # v is a new function
v(2,3);
                                    
 

For  Q  everything is similar.