@tomleslie In the text part of the help on plots:-arrow, the option  color  is not really mentioned, but in the examples it is used. Therefore, it is natural to think that it was simply forgotten to mention in the text part of the help. On the other hand, this option is also absent in the text part of  plots:-display  command.
If we follow your recommendation, then in order to plot even one arrow with the needed color, we should call plots:-display  command?...

@ssllys  I just noticed that  af_eq  is quadratic with respect to sigma. Solving this quadratic equation, we get 2 branches  curve1  and  curve2  (right and left). Here is an example of working with these branches  (I plot these branches in different colors):

af_eq:=(3*alpha[2](t[1], t[2])*A^3*(1/8)-A*sigma)^2+(1/4)*alpha[1](t[1], t[2])^2*A^2 = (1/4)*f0^2;
f0 :=0.1; alpha[1]:=0.0087;alpha[2]:=2.5871;
af_eq;
curve1,curve2:=solve(af_eq, sigma);
minimize(curve1, A=0.1..0.5, location);
A0:=rhs(%[2,1,1,1]);
with(plots):
P1:=implicitplot([sigma=curve2,sigma=curve1],sigma = -1.2 .. 2, A = 0 .. A0,  numpoints = 20000, axes = box, axesfont=[SYMBOL, 14],labels = [sigma, A], labelfont = [SYMBOL, 16],color=[red,green],tickmarks=[9,12],thickness=2):
P2:=implicitplot([sigma=curve2,sigma=curve1],sigma = -1.2 .. 2, A = A0 .. 2,  numpoints = 20000, axes = box, axesfont=[SYMBOL, 14],labels = [sigma, A], labelfont = [SYMBOL, 16],color=[red,blue],tickmarks=[9,12],thickness=2):
display(P1,P2);
                       

@Ronan 

fact := proc (expr::equation) options operator, arrow; frontend(gcd, [numer(lhs(expr)), numer(rhs(expr))])/frontend(gcd, [denom(lhs(expr)), denom(rhs(expr))]) end proc;
 
eq2 := (F+G)*R*sin(x)/sqrt(W-1) = (1/2)*M*R^2*tan(x)*(a/(R*sqrt(W-1)));
eq22:=convert(eq2, sincos);
 fact(eq22); eq22/fact(eq22);

@tdavid  If  rho(phi)  is known, then you must represent this function here. For this you can (after saving your document first) upload it here, using the bold green up-arrow in the MaplePrimes editor.

For example if a curve is an ellips  x^2/5^2+y^2/3^2=1  then there is no problems:

rho:=t->1.8/(1+0.8*cos(t));  # Polar equation of this ellips is  r=rho(t)  (the pole in the right focus)
L:=phi->Int(sqrt(rho(t)^2+diff(rho(t),t)^2), t=0..phi);
plot(L(phi), phi=0..2*Pi);
evalf(L(4));
evalf(L(5));
phi=fsolve(L(phi)=50);

@tzeng  Below is the corrected code:

restart;
Mphi := f->-I*(diff(f, phi));
Pr := f->-I*(diff(f, r));
rplus := f->r*exp(I*phi)*f;
rminus := f->r*exp(-I*phi)*f;
Pplus := f->exp(I*phi)*(Pr(f)+I*Mphi(f)/r);  
Pminus := f->exp(-I*phi)*(Pr(f)-I*Mphi(f)/r);
Rpp := f->(Pplus(f)+I*rplus(f));
Rmp := f->(Pplus(f)-I*rplus(f));
Rpm := f->(Pminus(f)+I*rminus(f));
mm := f->(Pminus(f)-I*rminus(f));
Psi00 := sqrt(1/Pi)*exp(-(1/2)*r^2);                        
Rpm(Rpp(Psi00));                
simplify(%);   
 

@Earl  Wonderful code! Thanks for this.

@Fereydoon_Shekofte  Here is a procedure for this.  N  is the number of the steps.

restart;
CubeWalk:=proc(N::posint)
local V1,V2,V3,V4,V5,V6,V7,V8,Faces,Colors,Cube,r,onward,backward,left,right,L,n,a,T,F1,F2,F3,F4,L1,L2,L3,L4,A; 
uses plottools,plots;

V1,V2,V3,V4,V5,V6,V7,V8:=[0,0,0],[0,1,0],[1,1,0],[1,0,0],[0,0,1],[0,1,1],[1,1,1],[1,0,1]:  # The vertices of the cube
Faces:=[[V1,V4,V8,V5],[V5,V6,V7,V8],[V2,V3,V7,V6],[V1,V2,V3,V4],[V3,V4,V8,V7],[V1,V2,V6,V5]]: # The list of the faces
Colors:=[green,red,magenta,blue,grey,gold]: # The list of the colors
Cube[0]:=display([seq(polygon(Faces[i],color=Colors[i]),i=1..6)]):

r:=rand(1..4);
onward,backward,left,right:=1,2,3,4;
L:=[V1,V2,V3,V4];
L1,L2,L3,L4:=[L[1],L[2]],[L[3],L[4]],[L[1],L[4]],[L[2],L[3]];

for n from 1 to N do
a:=r();

if a=onward then 
F1[n]:=t->rotate(Cube[n-1],t, L1);
Cube[n]:=rotate(Cube[n-1],-Pi/2, L1);
T:=textplot3d([0,0,1.5,typeset("Step ",n)], align='left', font=[times,16]);
A[n]:=animate(display@F1[n],[t], t=0..-Pi/2, background=T, paraminfo=false, frames=5); 
L:=map(p->p+[-1,0,0], L) else


if a=backward then 
F2[n]:=t->rotate(Cube[n-1],t, L2);
Cube[n]:=rotate(Cube[n-1],-Pi/2, L2);
T:=textplot3d([0,0,1.5,typeset("Step ",n)], align='left', font=[times,16]);
A[n]:=animate(display@F2[n],[t], t=0..-Pi/2, background=T, paraminfo=false, frames=5); 
L:=map(p->p+[1,0,0], L) else

if a=left then 
F3[n]:=t->rotate(Cube[n-1],t, L3);
Cube[n]:=rotate(Cube[n-1],Pi/2, L3);
T:=textplot3d([0,0,1.5,typeset("Step ",n)], align='left', font=[times,16]);
A[n]:=animate(display@F3[n],[t], t=0..Pi/2, background=T, paraminfo=false, frames=5); 
L:=map(p->p+[0,-1,0], L) else

if a=right then 
F4[n]:=t->rotate(Cube[n-1],t, L4);
Cube[n]:=rotate(Cube[n-1],-Pi/2, L4);
T:=textplot3d([0,0,1.5,typeset("Step ",n)], align='left', font=[times,16]);
A[n]:=animate(display@F4[n],[t], t=0..-Pi/2, background=T, paraminfo=false, frames=5);
L:=map(p->p+[0,1,0], L);

fi; fi; fi; fi;

L1,L2,L3,L4:=[L[1],L[2]],[L[3],L[4]],[L[1],L[4]],[L[2],L[3]]; 

od:

display([seq(A[k], k=1..N)], insequence, scaling=constrained, axes=normal);

end proc:

Example of use:

 CubeWalk(20);
                                    

Edit.    

@vv  Thanks for this workaround. Probably still need this replacement should be considered as a serious bug in Maple 2018.1 .

@one man 

1. In your document   cube_without_slipping.mw , the cube moves only in one direction (from left to right), in contrast to the animation presented at MaplePrimes.

2. Below is another code that does the same animation. Using  plottools:-rotate  command allows you to significantly simplify the code. Using  plots:-animate  instead  of  plots:-display  with  insequence  is also more convenient because it automates the creation of individual frames. Also programmatically made the  movement to the left.

3. For some unknown reason, the code below does not work in Maple 2018.1, but works in Maple 2015 and Maple 2017

restart; 
with(plottools): with(plots):
V1,V2,V3,V4,V5,V6,V7,V8:=[0,-1,0],[0,0,0],[1,0,0],[1,-1,0],[0,-1,1],[0,0,1],[1,0,1],[1,-1,1]:  # The vertices of the cube
Faces:=[[V1,V4,V8,V5],[V5,V6,V7,V8],[V2,V3,V7,V6],[V1,V2,V3,V4],[V3,V4,V8,V7],[V1,V2,V6,V5]]: # The list of the faces
Colors:=[green, red,RGB(1, 0, 4),blue,grey,gold]: # The list of the colors
Cube[0]:=display([seq(polygon(Faces[i],color=Colors[i]),i=1..6)]):

for n from 1 to 7 do
F[n]:=t->rotate(Cube[n-1],t, [[0,n-1,0],[1,n-1,0]]):
Cube[n]:=rotate(Cube[n-1],-Pi/2, [[0,n-1,0],[1,n-1,0]]):
A[n]:=animate(display,[F[n](t)], t=0..-Pi/2, paraminfo=false);
od:

for m from 6 to 0 by -1 do
G[m]:=t->rotate(Cube[m+1],t, [[0,m,0],[1,m,0]]):
B[m]:=animate(display,[G[m](t)], t=0..Pi/2, paraminfo=false);
od:

C1:=display([seq(A[k], k=1..7)], insequence):
C2:=display([seq(B[k], k=6..0, -1)], insequence):
display([C1,C2], insequence, scaling=constrained, axes=normal);

@Ronan  In the MaplePrimes editor, click this button

@Ronan   See update to my answer.

@vv  Thank you for this example. It was quite unexpected for me. This is very instructive, because allows more critical attitude to  RootOf  representation for roots of equations.

@9009134  I do not know the answers to your questions.

@dumin 

f1:=8.044048-0.764286*x1-0.756746*x2+0.034524*x1*x1+0.022222*x2*x1+0.098413*x2*x2;
maximize(f1, x1=3..9, x2=1..7, location);

Addition.  If we consider  f1 on the whole plane  R^2, then the function is unbounded from above, i.e. can take any arbitrarily large values.

@Carl Love I think that the reason for the bug is related to the property of the function  Ei(x)  at  x=0. It would be better to add more  assuming epsilon>0  to my code.

Addition. I just noticed that the bug only occurs if the lower limit of the improper integral is 0. In other cases, the integral is calculated correctly, for example:

restart;
int(exp(-t)/(1-t), t = -0.1 .. infinity, CauchyPrincipalValue = true);
int(exp(-t)/(1-t), t = 0.1 .. infinity, CauchyPrincipalValue = true);
                                 0.7973339105
                                 0.5969990665