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);

Use functional assignment for the calculating  f   for different points.

See corrected file.

Download failforce_new.mw

Unfortunately, GraphTheory:-PrimsAlgorithm command returns a list of edges for the minimal spanning tree not in the same order as these edges are obtained by applying Prim's algorithm. Here is another procedure named  PrimsAlg  that returns a list of edges in the correct order. The procedure does not use any commands from GraphTheory package and works strictly according to the Prim's algorithm.The required parameter  L  is the list of the vertices of the graph given by the coordinates of points on the Euclidean plane. The optional parameter  N   is the number of the vertex with which we begin to build the minimal spanning tree ( by default N=1).

The code of the procedure:

PrimsAlg:=proc(L::listlist, N::posint:=1)
local n, A, dist, L1, V, U, i, p, E;
n:=nops(L);
A:=Matrix(n, (i,j)->evalf(sqrt((L[i,1]-L[j,1])^2+(L[i,2]-L[j,2])^2)));

dist:=proc(U::set(posint), V::set(posint))
local a, u, v, p;
a:=infinity;
  for u in U do
    for v in V do
      if A[u,v]<a then a:=A[u,v]; p:=[u,v,a] fi;
    od; 
  od;
p;
end proc:

L1:=convert(L, set):
V[1]:={N}: U[1]:= {$1..n} minus V[1];

  for i from 1 to n-1 do
    p:=dist(V[i],U[i]);
    V[i+1]:=V[i] union {p[2]}; U[i+1]:=U[i] minus {p[2]};
    E[i]:=p[1..2];
  od;

convert(E, list);

end proc:

Example of use for the list of 50 random points in the square  [0,30] x [0,30]:

L:=RandomTools:-Generate(listlist(integer(range = 0 .. 30), 53, 2)):  # Generating 53 random points
M:=convert(convert(L,set)[1..50], list):  #  Selection of 50 different points
M1:=subs({seq(n=M[n], n=1..50)}, PrimsAlg(M)):  # Consecutive list of added edges

M1 := [[[0, 3], [2, 1]], [[0, 3], [2, 7]], [[2, 7], [2, 8]], [[2, 8], [0, 11]], [[2, 8], [4, 11]], [[4, 11], [6, 8]], [[4, 11], [8, 12]], [[8, 12], [9, 16]], [[9, 16], [8, 17]], [[8, 17], [8, 20]], [[8, 20], [8, 21]], [[8, 21], [10, 24]], [[10, 24], [8, 26]], [[8, 26], [8, 29]], [[8, 12], [12, 11]], [[12, 11], [13, 12]], [[8, 17], [4, 16]], [[4, 16], [1, 17]], [[1, 17], [2, 21]], [[8, 26], [4, 25]], [[4, 25], [4, 28]], [[10, 24], [14, 23]], [[14, 23], [14, 22]], [[14, 23], [17, 23]], [[17, 23], [16, 27]], [[16, 27], [15, 28]], [[15, 28], [13, 30]], [[16, 27], [19, 29]], [[17, 23], [21, 22]], [[12, 11], [16, 8]], [[16, 8], [17, 8]], [[17, 8], [18, 8]], [[18, 8], [20, 9]], [[20, 9], [20, 13]], [[20, 13], [22, 13]], [[22, 13], [23, 16]], [[18, 8], [20, 4]], [[20, 4], [18, 2]], [[18, 2], [16, 1]], [[20, 4], [23, 4]], [[23, 4], [27, 4]], [[27, 4], [27, 1]], [[27, 1], [29, 0]], [[27, 1], [24, 0]], [[27, 4], [30, 7]], [[22, 13], [27, 13]], [[23, 16], [28, 19]], [[28, 19], [29, 19]], [[19, 29], [27, 28]]]

Visualization:

A:=plot(M, style=point, color=blue, symbol=solidcircle, symbolsize=15):
B:=seq(plot(M1, color=red, thickness=2), i=1..49):
plots:-display(A, B);  # The minimal spanning tree for M

Addition - animation of the building of the minimal spanning tree  (the code below does not depend on the previous one):

M1 := [[[0, 3], [0, 3]], [[0, 3], [2, 1]], [[0, 3], [2, 7]], [[2, 7], [2, 8]], [[2, 8], [0, 11]], [[2, 8], [4, 11]], [[4, 11], [6, 8]], [[4, 11], [8, 12]], [[8, 12], [9, 16]], [[9, 16], [8, 17]], [[8, 17], [8, 20]], [[8, 20], [8, 21]], [[8, 21], [10, 24]], [[10, 24], [8, 26]], [[8, 26], [8, 29]], [[8, 12], [12, 11]], [[12, 11], [13, 12]], [[8, 17], [4, 16]], [[4, 16], [1, 17]], [[1, 17], [2, 21]], [[8, 26], [4, 25]], [[4, 25], [4, 28]], [[10, 24], [14, 23]], [[14, 23], [14, 22]], [[14, 23], [17, 23]], [[17, 23], [16, 27]], [[16, 27], [15, 28]], [[15, 28], [13, 30]], [[16, 27], [19, 29]], [[17, 23], [21, 22]], [[12, 11], [16, 8]], [[16, 8], [17, 8]], [[17, 8], [18, 8]], [[18, 8], [20, 9]], [[20, 9], [20, 13]], [[20, 13], [22, 13]], [[22, 13], [23, 16]], [[18, 8], [20, 4]], [[20, 4], [18, 2]], [[18, 2], [16, 1]], [[20, 4], [23, 4]], [[23, 4], [27, 4]], [[27, 4], [27, 1]], [[27, 1], [29, 0]], [[27, 1], [24, 0]], [[27, 4], [30, 7]], [[22, 13], [27, 13]], [[23, 16], [28, 19]], [[28, 19], [29, 19]], [[19, 29], [27, 28]]]:
M:={op~(M1)[ ]}:
A:=plot(M, style=point, color=blue, symbol=solidcircle, symbolsize=15):
F:=p->plots:-display(seq(plottools:-line(M1[i][], color=red, thickness=2), i=1..p));
plots:-animate(plots:-display,['F'(p)], p=1..50, background=A, frames=51);

Prim.mw

Edit.

You can use  textplot  command instead of  legend  option to place the legends in any desired place. Of course, this will require a few more efforts.

Example:

with(plots):
A:=plot([sin(x), cos(x)], x=-Pi..3*Pi, color=[red,blue]):
B:=textplot([[3.7,0.8,sin(x)], [-1.8,0.5,cos(x)]], font=[times,bold,
14]):
C:=textplot([[3,1,__, color=red], [-2.5,0.7,__,color=blue]], font=[times,bold,24]):
display(A, B, C, size=[950,200], scaling=constrained);

Use  PolyhedralSets  package. See the help on this package for details. Here is an example of the plotting of your first set (all inequalities must be non-strict):

with(PolyhedralSets):
P:=PolyhedralSet({40 <= y, x <= y-40, z <= y-40});
Plot(P, color=khaki, axes=normal, view=[-30..60, 0..90, -40..50], orientation=[40,70]); 

We see that this is an unbounded trihedral body.
 

A lot of syntax errors. Also note that  implicitplot  command only works for two variables. To examine the dependence of the plot on the third variable (I took  p  as the third variable), you can use  animate  or Explore  commands:

restart; 
with(plots):
omega[0]:=1: epsilon[1]:=1.04493: epsilon[2]:=0.93259: delta:=0.02:
eq:=(omega[0]*a-a*Omega^2+(3/4)*epsilon[2]*a^3)^2+(delta*Omega*a)^2 = p^2;
animate(implicitplot, [eq, Omega = 0.5 .. 2, a = 0 .. 5, color = blue, thickness = 2, gridrefine = 3], p=0..0.2, frames=120);

# Or
Explore(implicitplot(eq, Omega = 0.5 .. 2, a = 0 .. 5, color = blue, thickness = 2, gridrefine = 3), parameters=[p=0...1.], numframes =120);

Edit.