F := k -> sum(i^k, i=1..n);

seems to be faster.

Check this FAQ:

https://desktopfaq.maplesoft.com/customer/en/portal/articles/2840902-maple-does-not-start-or-takes-a-long-time-to-load-

See:

https://www.mapleprimes.com/posts/207910-A-Compact-Sudoku-Solver-Via-Groebner

Probably the designers considered that the color, linestyle, thickness are enough for the legend (which I also do).

Anyway, the symbols seem to be useless here because their exact position is not known.

You have the Physics package loaded. It redefines many things, including diff.
The conversion works as expected if you don't use the package (which is not needed in the worksheet).

If you have the curve in polar coordinates, then DirectSearch or fsolve is not needed.
Something like this (for a modified cardioid):

r:=t->2*(1-cos(t))+1:
plots:-animate(plot, 
[[[r(t-u)*cos(t),r(t-u)*sin(t),t=0..2*Pi] ,
[[0,r(Pi/2-u)*sin(Pi/2)],[0,2+r(Pi/2-u)*sin(Pi/2)]]   ],thickness=[1,6]],
u=0..2*Pi, scaling=constrained);

You cannot use plot3d(f(x,y)=0.1, x=0..1,y=0..1);
Probably you want plots[implicitplot](f(x,y)=0.1, x=0..1,y=0..1);

For z=0..1  just use the view option.
Compare:

plot3d(x^2+y^2, x=-1..1, y=-1..1);
plot3d(x^2+y^2, x=-1..1, y=-1..1, view=0..1);

In my opinion the student should be taught to use a single simple line of code

s := n -> s(n-1)^2 - 2;   s(1) := 4;

defining his own sequence:
s(4);
      37634
 seq(s(n), n=1..5);  
      4, 14, 194, 37634, 1416317954

See a full explanation here: https://en.wikipedia.org/wiki/Zero_to_the_power_of_zero

evalindets treats the type `*` using the internal representation rather than using op.
For example 5*x^3*y  and a*x^3*y  have distinct internal representations (check with dismantle).
The next procedure has not this problem.

M:=proc(f::polynom,x,y)
local i,t, a:=coeffs(f,[x,y],t);
add(a[i]* m[degree(t[i],x),degree(t[i],y)], i=1..nops([a]));
end;

M(5*x^3*y+13*a*y^2+z, x,y);

You have a mess in your arguments. You substitute f(t,s)  but keep  s,x?
Actually it is not clear which is the intended  function  subject to invhilbert.

Sierpinski:=proc(n,cir:=true)
local L:=[[0,sqrt(3)],[-1,0],[1,0],[0,sqrt(3)]], S;
S:=proc(a,b,c)
  L:=L,[(b+c)/2,(a+c)/2,(a+b)/2,(b+c)/2];
  if (c[1]-b[1])>2^(-n+2) then
     S(a,(a+b)/2,(a+c)/2);
     S((a+b)/2,b,(c+b)/2);
     S((a+c)/2,(b+c)/2,c) fi;
  end;
S([0,sqrt(3)],[-1,0],[1,0]);
plot([L],axes=none,scaling=constrained,thickness=0,title="Sierpinski "||(nops([L])),color=blue);
plots:-display(%,`if`(cir,plot([L], style=point,symbol=circle,symbolsize=10,color=blue),NULL));
end:
Sierpinski(4); Sierpinski(7,false);

    a = b mod n;

is interpreted as

   a = (b mod n)
or (by default) as
    a = modp(b, n)

see, the help page ?operators,precedence

So, 

3 = 0 mod 3;
                             3 = 0
0 = 3 mod 3;
                             0 = 0