with(plots):with(plottools):
p1:=plot( x^3,x=-1..1,thickness=20,color=red):
p2:=plot(-x^3,x=-1..1,thickness=20,color=blue):
p3:=display(rectangle([0.5, 1],[0.75,-1],color=green)):
p4:=display(rectangle([-1,-0.1],[1,0.1],color=yellow)):

display(p1,p2,p3,p4); # order = 4312

#  display(p3,p4,p1,p2); # order = 4312 (the same)

It seems that the rectangles are plotted first, in reversed order, and then the curves, in direct order.
Has someone an explanation?

There are many conventions for the Euler angles or other angles used to define o rotation of a 3d plot.
In Maple these angles are in the plot option orientation, but I think that the help page is not correct about them.
The same info appears in a worksheet (see ?rotateplot), so I am even more intrigued. [Note that many authors also switch phi and theta in spherical coordinates].

The help file says:

orientation=[theta, phi, psi]
This orientation specified by these angles is obtained by rotating the plot
1. psi about the x-axis,
2. then phi about the (transformed) z-axis, and
3. then theta about the (transformed) y-axis.
  These angles, given in degrees, are the Euler angles for the transformation matrix, using the axes specified. The angle psi is optional and is assumed to be 0 if not given. If the orientation option is not specified, the default value used is [55, 75, 0].

After some tests it seems that y and z should be switched, i.e. keeping the names (and order) for the angles ==>
1. psi about the x-axis,
2. then phi about the (transformed) y-axis, and
3. then theta about the (transformed) z-axis.
 
Am I right?

It seems to be a bug here:

restart;
with(Iterator):
M := CartesianProduct([1,2],[a,b],[c,d,e], rank=3):
n:=0:  for v in M do n:=n+1: print(n,v); end do:

              1, [2   b   c]
              2, [2   b   d]
Error, (in unknown) improper op or subscript selector

# for rank=9  ==>  [list, b, e]     ??

g:= a -> int(f(x+a),x=a..2*a):
eval(g(x),x=1);
                    int(f(2 x), x = 1 .. 2)
eval(g(z),z=1);
                   int(f(x + 1), x = 1 .. 2)

eval is advertised as smart, but it's not enough!


(Related to a recent answer of Carl Love about bound variables.)

Edit: shorter version

g:= a -> int(sin(sin(x+a)),x=a..2*a):
evalf(eval(g(x),x=1))=evalf(eval(g(z),z=1));
                  0.1052070507 = 0.5294405453

minimize((x+y)^2+cos(y)^2, x=-4..4, y=-2..2);  # cos(2)^2    should be 0
minimize((x+y)^2+cos(y), x=-4..4, y=-4..4);      
   Error, (in unknown) mapped procedure in `ormap' must return true or false
minimize((x+y)^2+cos(y));                                   # infinity   ?!  should be -1

Edited. Corrected second example