@Michael_Watson 

@tomleslie Maybe Ronan wants some global variables to be used enywhere. Of course it's not the best practice but it's simple and useful sometimes.

@Detlef Hoyer It would break some standard constructs.
E.g. (exp - 1)(x)  will not be any longer the expected  exp(x) - 1.

@mmcdara The procedure returns true if the two varieties are the same. For your second example it gives correctly "false" because the first variety is a line and the second is a plane.
Note that it works for more general varieties, e.g.
EQSYS( {x+y+z=1}, {x=u+1, y=v+1, z=-u-v-1}, param={u,v});   # true

The procedure returns an error if one of the varieties is empty, but it is easy to change this:

EQSYS:=proc(S1::set(`=`), S2::set(`=`), {vars:={x,y,z}, param:={t}}) 
  local A:=eliminate(S1,param)[2], B:=eliminate(S2,param)[2],
        SA:=solve(A,vars), SB:=solve(B,vars);
  if SA=NULL then SA:={} fi; if SB=NULL then SB:={} fi;
  evalb( eval(B,SA) union eval(A,SB) = {0});
end:

@mmcdara For b) the invariant is  L(M(a,b), N(a,b)) = L(a,b).
So, L(a[n], b[n]) = L(a[0], b[0])   and let n --> oo.

@GunnerMunk You should not use the document mode for teaching. It is used mainly for presentations.

@dharr I don't think there are many such bugs.
The series structure is special and usually it is not part of an algebraic expression; actually the arithmetic of series works only inside series. For the second example, the user just has to know about the atomic expressions.
So, one may use (in general) the type system without knowing all the details of Maple internals (e.g. DAGs).

@Joe Riel  Because

type(x^2, identical(x)^2);  and   type(x^2, identical(x^2));

are both true.

You have only two variables (eta and u). So, what kind of 3D do you mean?

The legends are separated from the plot area in Maple. You may use textplot etc if you really want this.

@mmcdara Yes, I also use this often. Vote up.
Unfortunately it may fail even for simple expressions.

J:=Int(Heaviside(1-x^2-y^2)*x^2, [x=-1..1, y=-1..1]):
K:=Int(Heaviside(1-x^2-y^2)*x^2, [x=-infinity..infinity, y=-infinity..infinity]):
value(J) = evalf(J);
#                        0 = 0.7853981634
identify(%);
#                        0 = Pi/4
value(K);
#                          undefined
evalf(K);
#                          0.7853981634

Note that even int(Heaviside(1-x^2-y^2), [x=-1..1, y=-1..1])  fails (it gives 0 too).

@mmcdara It was approx 70 sec,  almost all for the second eliminate (for I2).

@mmcdara  If `union`(S) = {a1,...,an}  and  nops(S) >> n,  probably a good candidate for F(S) is  { {}, {a1}, ..., {an} }.

@tomleslie I don't see other solution if the present functionality is maintained. The default labels could be `1`, `2`,... or maybe 1., 2., ... (non-integer). Or, it may remain as it is, with a warning.

@Earl If the contact points are allowed to be anywere on the lines, any conic is possible.

@Carl Love I ment f (z,A,phi) has values arbitraryly close to 0 in any nbd of complex infinity (complex infinity being an essential singularity of f(., A, phi)). So, it will be impossible to find numerically all the roots, if a bounded region is not known.