Add the following line to your code:
map(simplify~, %[2], zero);

See corrected file:

ANALYTIC_1_1.mw

Do this in polar coordinates:

plot3d(eval([x, y, x^3-2*x^2*y], [x=r*cos(phi), y=r*sin(phi)]), r=0..1, phi=0..2*Pi, axes=normal, numpoints=5000);

In Cartesian coordinates, this can also be done, but the quality of plotting is slightly worse:

plot3d(x^3-2*x^2*y, x=-1..1, y=-sqrt(1-x^2)..sqrt(1-x^2), axes=normal, numpoints=5000);

Edit.
 

You have come up with an intricate and very inefficient method of solving the problem, in which you made a lot of mistakes. Such tasks are instantly solved using seq command:

make_PSI:=nplanets->[seq(2*Pi*k/nplanets, k=0..nplanets-1)]:

Example of use:

make_PSI(4);

                             [0, (1/2)*Pi, Pi, 3*Pi*(1/2)]
 

with(plots):
Sys:=[x(t) = cos(t) + t*sin(t), y(t) = sin(t) - t*cos(t)]:
C:=spacecurve([rhs~(Sys)[ ], 0], t=0..Pi/2, axes=normal, color=red, thickness=4): # The curve in 3d
S:=plot3d(eval([x(t), y(t)*cos(phi), y(t)*sin(phi)], Sys), phi=0..2*Pi, t=0..Pi/2): # The surface of rotation
display(C, S, style=surface, scaling=constrained, labels=[x,y,z], orientation=[70,70]);  # Alltogether

Addition - animation of the rotation:

RC:=phi->plottools:-rotate(C, phi, [[0,0,0],[1,0,0]]):
RS:=phi->plot3d(eval([x(t),y(t)*cos(s),y(t)*sin(s)],Sys), s=0..phi, t=0..Pi/2, style=surface):
animate(display, ['RC'(phi), 'RS'(phi)], phi=0..2*Pi, scaling=constrained, orientation=[70,70], lightmodel=light1, frames=90);

SurfaceOfRotation.mw

eq := diff(x(t), t) = x(t) + 1:
indets(eq, {name, function(name)});

                                {t, x(t)}

Maybe this is a bug (the first version with the loop)?

Here is a workaround:

nonIdMaps := [ ]:
F:=[x -> x,x -> 2*x,x -> 3*x]:
for i from 1 to 3 do
F[i], F[i](y);
if F[i](y) <> y then nonIdMaps := [nonIdMaps[],F[i]] end if;
end do;
nonIdMaps, map(f -> f(y),nonIdMaps);

Try

labels = [x, `#mover(mi("sigma"),mo("&circ;"))`[y]]
 

# Or

labels = [x, conjugate(sigma[y])]

Edit.

See help on  mtaylor  command.

Example. Should be an initial condition for the numerical solution:

sol:=dsolve({diff(y(x),x) = -2*x*y(x) + 1, y(0)=2}, numeric, method=classical, stepsize=0.1); # The numerical solution in the form of a procedure 
plots:-odeplot(sol, [x,y(x)], x=0..1); # The plotting of the solution in the specific range 

# Finding the values of the function y(x) at individual points
eval(y(x), sol(0.1));
eval(y(x), sol(0.15));
eval(y(x), sol(0.2));

We see that to find the values of the function at intermediate points Maple uses linear interpolation in this method.

dsolve.mw

Edit.

See  desired_form_in_output_new.mw

Example:

Points:=[[1,1],[1,3],[2,2],[3,3],[3,1]]:  # The points in the form of a list of their coordinates
P:=plot(Points, style=point, color=red, symbol=solidcircle, symbolsize=20):  # The points
L:=plot(Points, color=blue, thickness=2):  # The line segments
plots:-display(P, L, view=[0..4, 0..4]);  # All together

Edit.

with(plots):

K:=9;
deG:=diff(theta(t),t,t) + mu*diff(theta(t),t)+K*sin(theta(t))= 0;
deL:=diff(theta(t),t,t) + mu*diff(theta(t),t)+K*theta(t)= 0;
Iv:=theta(0)=0.75, D(theta)(0)=2.0;
dom1:=t=0..10;
soln1a:=dsolve({eval(deL,mu=0),Iv});

soln1b_1:=dsolve({eval(deL,mu=1),Iv});

gr1c:=multiple(plot,[rhs(soln1a),dom1,color=blue],[rhs(soln1b_1),dom1,color=purple]);
 

Addition.

In fact, there is no need to call plots package and use multiple command. It's simpler and smarter to use the standard  plot command, and write expressions to be plotted as a list:

restart;
K:=9;
deG:=diff(theta(t),t,t) + mu*diff(theta(t),t)+K*sin(theta(t))= 0;
deL:=diff(theta(t),t,t) + mu*diff(theta(t),t)+K*theta(t)= 0;
Iv:=theta(0)=0.75, D(theta)(0)=2.0;
dom1:=t=0..10;
soln1a:=dsolve({eval(deL,mu=0),Iv});
soln1b_1:=dsolve({eval(deL,mu=1),Iv});

plot([rhs(soln1a),rhs(soln1b_1)], dom1, color=[blue,purple]);  # Plotting by the usual  plot  command

Edit.
 

Some commands from  Student[VectorCalculus]  subpackage solve a number of problems of classical differential geometry of curves, for example  Binormal ,  PrincipalNormal ,  Curvature ,  RadiusOfCurvature,  Torsion  and some others.

 
 ,

plot1:=plot(x, x=0..3, 0..7, color=red):
plot2:=plot(x^2, x=0..3, 0..7, color=blue):
plot3:=plot(x^3, x=0..3, 0..7, color=green):

plots:-animate(plots:-display,[[plot([[0,0]]), plot1$20, plot2$20, plot3$20][1..ceil(n)], scaling=constrained], n=1..60);

Addition.

Here is another way in which the technique from this post is used:

with(plots):
plot1:=animate(plot,[x, x=0..a, 0..7, color=red, thickness=2], a=0..3, frames=30):
plot2:=animate(plot,[x^2, x=0..a, 0..7, color=blue, thickness=2], a=0..3, frames=30):
plot3:=animate(plot,[x^3, x=0..a, 0..7, color=green, thickness=2], a=0..3, frames=30):
display([plot1, display(op([1,-1,1],plot1),plot2), display(op([1,-1,1],plot1),op([1,-1,1],plot2),plot3)], insequence, scaling=constrained);