Unfortunately, the developers did not provide for legends such as you want, so in the example below I made legends using the packages  plottools  and  plots :

restart;
with(plottools):  with(plots):
A:=plot(sin(x), x=0..2*3.14, style=pointline, symbol= diamond, color=red, symbolsize=11, numpoints=40, adaptive=false):
B:=plot(cos(x), x=0..2*3.14, style=pointline, symbol= solidcircle, color=blue, symbolsize=11, numpoints=40, adaptive=false):
L1, L2:=line([2.75,0.85],[3.2,0.85],color=red), line([3.65,0.85],[4.1,0.85],color=blue):
P1:=plot([seq([2.75+i*(3.2-2.75)/3,0.85], i=1..2)],style=point,symbol= diamond, color=red, symbolsize=11):
P2:=plot([seq([3.65+i*(3.2-2.75)/3,0.85], i=1..2)],style=point,symbol= solidcircle, color=blue, symbolsize=11):
T:=textplot([[3.4,0.85,sin(x)],[4.3,0.85,cos(x)]]):
P:=polygon([[2.7,0.7],[2.7,1],[4.5,1],[4.5,0.7]], style=line, thickness=0):
plots:-display(A, B, P, L1, L2, T, P1, P2, scaling=constrained, size=[800,300]);

Addition. This is done in Maple 2017.3. In earlier versions, you may have to specify graphs for  style=line  and  style=point  separately.

restart;
L:=[[0,0,0], [1,0,0], [1,1,0], [0,1,0], [0,0,1], [1,0,1], [1,1,1], [0,1,1]]:
T:=combinat:-choose(L, 4):
Dist:=(p,q)->sqrt((p[1]-q[1])^2+(p[2]-q[2])^2+(p[3]-q[3])^2):
P:=v->combinat:-choose(v, 2):
seq(`if`(nops(convert((Dist@op)~(P(t)),set))=1, t, NULL), t=T);  # Final result (two solutions)

             [[0, 0, 0], [0, 1, 1], [1, 0, 1], [1, 1, 0]], [[0, 0, 1], [0, 1, 0], [1, 0, 0], [1, 1, 1]] 

A visualization:

restart;
with(plottools): with(plots): with(combinat):
A:=cuboid([0, 0, 0], [1, 1, 1], style=line, color=black):
L:=[[[0, 0, 0], [0, 1, 1], [1, 0, 1], [1, 1, 0]], [[0, 0, 1], [0, 1, 0], [1, 0, 0], [1, 1, 1]]]:
P1, P2:=choose(L[1], 2), choose(L[2], 2):
B:=seq(line(op(P1[i]), color=red, thickness=3), i=1..6):
C:=seq(line(op(P2[i]), color=blue, thickness=3), i=1..6):
display(A, B, C, axes=none);

Edit.

x:=(t,vter)->vter*(1-exp(-t/vter)):
y:=(t,vter)->(1+vter)*vter*(1-exp(-t/vter))-vter*t:
plot([seq([x(t,v),y(t,v),t=0..2.5], v=[0.3,1,3]), [limit~([x(t,v),y(t,v)], v=infinity)[], t=0..2.5]], color=[brown,red,black,blue], thickness=2, view=[0..2.5,-0.3..0.5], numpoints=10000, legend=(vter=~[0.3,1,3,infinity]), legendstyle=[font=[times,14],location=right], size=[600,400]);
 

Edit.

Use  minimize  command instead of extrema. See

2_new.mw

As an alternative you can use  Optimization:-Minimize  command:

Optimization:-Minimize(V, x = 0 .. 1); 

Do  restart  first and everything will work properly:

restart; 
solve(a = b+5, b); 
solve(a*x^2+b*x+c = 0, x);
                                                  

Addition. Directly in your document I get a very strange result:

I think that you should create a new document that is free of any presetting.

It is true for any real  a  and  b  if  b>0 :

is(arctan(a,b)=arctan(a/b)) assuming real, b>0;
                                               true

Download q1nwtnnonlinearsys_new.mw

Function  procedure solves the problem. Required parameters of the procedure: F is an expression of the variable  x  or a procedure, R is a list of the end points of a segment (the domain of the function), N is the number of points (the ends of the segment  R  are included in this number so  N>=2).  Optional parameters: the equality  output=list  (by default) or   output=plot , Opts  is a list (additional options for plotting). 

Function:=proc(F::{algebraic,procedure}, R::list, N::posint, Eq::`=`:=(output=list), Opts::list:=[])
local f, a, b, h, L;
if F::algebraic then f:=x->F else f:=F fi;
a,b:=op(R); h:=(b-a)/(N-1);
L:=[seq(evalf[3]([a+i*h,eval(f(x),x=a+i*h)]), i=0..N-1)];
print(L);
if Eq=(output=plot) then 
plot([f(x), L], x=a..b, color=[red,blue], op(Opts)) fi;
end proc:

Example of use:

Function(x^2, [-2,2], 9, output=plot, [scaling=constrained, size=[500,500]]);

You can use  plots:-contourplot  command
See also this post  https://www.mapleprimes.com/posts/202222-Contour-Curves-With-Labels

Examples of use:

restart:
f(x,y):=exp(-(x^2+y^2)/3)*cos(2*x*y):
plots:-contourplot(f(x,y), x=-2.5..2.5, y=-2.5..2.5, contours=[seq(-0.4 .. 1, 0.1), 0.95] ,filledregions=true, coloring=["White","DarkViolet"], numpoints=50000); 

ContoursWithLabels(f(x,y), -2.5 .. 2.5, -2.5 .. 2.5, {seq(-0.4 .. 1, 0.1), 0.95},Coloring = [colorstyle = HUE, colorscheme = ["White","DarkViolet"], style = surface]);

  Edit.      

I tried to open the attached file by selecting  Maple Input , but my computer freezes after that.
Try to restore your worksheet as follows:

Menu  File -> Recent Documents -> Restore Backup

plot([BesselJ(0,x), BesselJ(1,x)], x=0..10, color=[red,blue]);

See help on  ?Bessel
 

In order to insert a text into each frame for this animation method, we need to create a new sequence of frames using  plots:-display  command:

restart;
with(plots):
N := 92: 
A := [seq(plot([[.85*sin(t)^3-2+1.25*i/N, .85*(13*cos(t)*(1/15)-(1/3)*cos(2*t)-2*cos(3*t)*(1/15)-(1/15)*cos(4*t)), t = 0 .. Pi*i/N], [-.85*sin(t)^3-2+1.25*i/N, .85*(13*cos(t)*(1/15)-(1/3)*cos(2*t)-2*cos(3*t)*(1/15)-(1/15)*cos(4*t)), t = 0 .. Pi*i/N], [sin(t)^3+2-1.25*i/N, 13*cos(t)*(1/15)-(1/3)*cos(2*t)-2*cos(3*t)*(1/15)-(1/15)*cos(4*t), t = 0 .. Pi*i/N], [-sin(t)^3+2-1.25*i/N, 13*cos(t)*(1/15)-(1/3)*cos(2*t)-2*cos(3*t)*(1/15)-(1/15)*cos(4*t), t = 0 .. Pi*i/N]], color = red, thickness = 5, view = [-2.3 .. 2.4, -1.4 .. 1]), i = 1 .. N)]:
 
A1:=[seq(display(textplot([[-1.9+(i-1)*1.2/N, 0,`Betty`],[1.9-(i-1)*1.2/N, 0,`Alf`]], font=[TIMES,16]), A[i]), i=1..N)]: 
 
display(A1, insequence = true, scaling = constrained, axes=none, title ="Alf & Betty", titlefont=[TIMES,18] );

verify(2*log(x), log(x^2)) assuming x>0;                          
evalb(2*log(x) = log(x^2)) assuming x>0;

                                             true
                                             true

restart:
h1:=z->1-(delta2/2)*(1 + cos(2*(Pi/L1)*(z - d1 - L1))):
h2:=z->1-(delta2/2)*(1 + cos(2*(Pi/L2)*(z - d2 - L2))):
h3:=z->1+(delta2/2):
K1:=((4/h1(z)^4)-(sin(alpha)/F)-h1(z)^2+Nb*h1(z)^4):
K2:=((4/h2(z)^4)-(sin(alpha)/F)-h2(z)^2+Nb*h2(z)^4):
K3:=((4/h3(z)^4)+(cos(alpha)/F)+h3(z)^2+Nb*h3(z)^4):
lambda1:=Int(K1,z=0..0.2):
lambda2:=Int(K2,z=0.2..0.4):
lambda3:=Int(K3,z=0.4..0.6):
lambda:=lambda1+lambda2+lambda3:
 
F:=0.3:
L1:=0.2:
d1:=0.2:
d2:=0.2:
L2:=0.3:
alpha:=Pi/6:
plot( [seq(eval(lambda, Nb=j), j in [0.1,0.2,0.3])], delta2=0.02..0.1);

Try adding this condition before  od;
if  not type(prO[i], float) then break fi;