restart;
S := {0.2, 0.5}:
R, Q := op(S);

R;
Q;

# Or

restart;
S := {0.2, 0.5}:
R:=S[1];
Q:=S[2];
R, Q;
 

Should be  de1 :=  ...;  de2 :=  ... ;   instead of    de1 =  ...;  de2 =  ... ;

Miquel  procedure allows you to visualize Miquel's theorem about five circles. The user himself arbitrarily sets the configuration parameters (any 5 numbers in the range  0..2*Pi). These numbers are the polar angles of the intersection points of blue circles with a stationary red circle (radius 1). The centers of all the blue circles lie on the red circle. Of course, this does not prove the theorem in the strict sense of the word, but it allows everyone to see clearly its fairness. The procedure does not use any special geometric packages, only  plottools  package for specifying lines and circles and  plots:-display command for plotting everything.

Download 111.mw

Your objective function  z  depends only on variables  x[i,j] , and the set of constraints also depends additionally on y[i]  and  u[i]  variables, which actually play the role of parameters (and the result should depend on these parameters). I think that Maple simply does not have the tools for solving similar problems (I do not know if such tools exist at all). To get a solution you must specify the values of these parameters.

restart;
eu := -2*n^2*B^2+2*n*B^2*p-2*n*B^2*p^2+n*B^2-B^2*p+B^2*p^3+n*c*d;
Eq1:=u=`p^2`+n-p:  Eq2:=v=`p^3`+n-p:
expand(subs([p^2=solve(Eq1,`p^2`), p^3=solve(Eq2,`p^3`)], eu)):
collect(%, [u,v]);
subs([subs(`p^2`=p^2, Eq1), subs(`p^3`=p^3, Eq2)], %);

Try this option:

restart;
p:=(a+b+c)/2: S:=sqrt(p*(p-a)*(p-b)*(p-c)):
Sol:=(hA,hB,hC)->allvalues(solve({hA=2*S/a, hB=2*S/b, hC=2*S/c,a>0,b>0,c>0}, {a,b,c}));

Example of use:

Sol(sqrt(3)/2, sqrt(3)/2, sqrt(3)/2);  
                                                                {a = 1, b = 1, c = 1}

I do not know whether this will work in Maple 7, because I do not have it.

Animation with legends:

f4 := seq(display(T || i, axes = boxed, size = [400, 320], titlefont = ["Arial", 12, Bold], labelfont = ["Arial", 13, Bold], labeldirections = [horizontal, vertical]), i = 1 .. nops(HA)): 
animate(display, [f4[round(i)], textplot([0.7, 0.8, cat(`Plot  `, `of  `, T, round(i))], font = [times, bold, 14])], i = 1 .. 4, frames = 60, paraminfo = false);

# Or

L := [0, 0.5, 1, 1.5]:
f4 := seq(display(T || i, axes = boxed, size = [400, 320], titlefont = ["Arial", 12, Bold], labelfont = ["Arial", 13, Bold], labeldirections = [horizontal, vertical]), i = 1 .. nops(HA)): 
animate(display, [f4[round(i)], textplot([0.7, 0.8, delta = L[round(i)]], font = [times, bold, 14])], i = 1 .. 4, frames = 60, paraminfo = false);

Edit.  

Expr:=(M/(L*T))^a*(M/L^3)^b*(L^2)^c;
simplify(Expr, symbolic);  # The simplest way
simplify(Expr) assuming M>0, L>0, T>0;  # More reliable way

I did not understand the code you wrote and therefore wrote a new code. I suggested that you want to animate the rotation of the ball around some axis specified by points  x  and  y . The coordinates of these points in the form of two lists are the parameters of  Rot  procedure:

Rot:=proc(x::list, y::list)
local Sph;
uses plots, plottools;
Sph:=plot3d(1, t=0..2*Pi, s=0..Pi, coords=spherical);
animate(display, [rotate(Sph, alpha, [x,y]), axes=normal, view=[(-1.35..1.35)$3]], alpha=0..2*Pi, frames=90, background=display(line(x,y), linestyle=3, color=red, thickness=2));
end proc:

Example of use:

Rot([0,0,0], [1,2,2]);

                                         
PS.   Of course, instead of a ball, you can take any other object to rotate it around the axis specified by points  x ,  y .  To do this, simply make an obvious change to the procedure code. You can also remove the rotated object from the body of the procedure and make it as the parameter of the procedure. Then you do not need to change the procedure code at all.
 

Edit.                     

Specify the values of the parameters and calculate numerically as in the example below. The imaginary part can be dropped:

Download Karim_new.mw

I do not have Maple 7, but I have  Maple 2015 - 2017. But in them (in geom3d package) I did not find a simple way to specify an arbitrary polyhedron, except for a few special cases. Therefore, I propose to use plottools package, in which each face of the polyhedron is given by polygon command. Using  plots:-display command you get a polyhedron bounded by these faces. First, in this way I set the faces of the tetrahedron, at the base of which the regular triangle with side 1, the height I took equal to 3/2. Then using the Pythagorean theorem 2 times, I calculated the height of the pyramid with a square base (also with side 1), which has the same lateral edges as the tetrahedron. Then I set its faces in the same way. plots:-display command allows you to depict everything on a single plot.

with(plots): with(plottools): 
ht:=3/2:
Vt:=[[0,0,0], [1,0,0], [1/2,sqrt(3)/2, 0], [1/2, sqrt(3)/6, ht]]:  # The vertices of the tetrahedron
t1, t2, t3, t4:=polygon(Vt[1..3], color=brown), polygon(subsop(3=NULL,Vt), color=blue), polygon(subsop(2=NULL,Vt), color=yellow), polygon(Vt[2..4], color=green):  # All faces of the tetrahedron
st:=evalf(sqrt((sqrt(3)/3)^2+ht^2));  # Length of lateral edge of the tetrahedron
hs:=evalf(sqrt(st^2-(1/sqrt(2))^2));  # Height of a pyramid with the square base
Vs:=[[0,0,0], [-1,0,0], [-1,1,0], [0,1,0], [-1/2, 1/2, hs]]:  # The vertices of the square pyramid
s1, s2, s3, s4, s5:=polygon(Vs[1..4], color=brown), polygon(subsop(2=NULL,3=NULL,Vs), color=green), polygon(subsop(3=NULL,4=NULL,Vs), color=blue), polygon(subsop(1=NULL,4=NULL,Vs), color=yellow), polygon(subsop(1=NULL,2=NULL,Vs), color=red):  # # All faces of the square pyramid
display(t1, t2, t3, t4, s1, s2, s3, s4, s5, scaling=constrained, axes=normal, view=[-1.35..1.35,0..1.35,0..1.55], orientation=[70,75], lightmodel=light1);

                      


Edit. Fixed errors in the formulas for st and hs.

Your equation has an obvious conststant solution  y=ln(2) . If you want to get some other solution, you can solve it numerically with the specified initial conditions, as in the example below:

ode := diff(y(x), x, x)-exp(y(x))+2 = 0;
ics:=y(0)=0.6, D(y)(0)=0;
Sol:=dsolve({ode,ics}, numeric);
plots:-odeplot(Sol, [x,y(x)], x=0..5,  view=-10..2, color=red);

PS. A solution in the form of a series usually gives a good approximation only in sufficient proximity to the starting point. Compare the two solutions:

ode := diff(y(x), x, x)-exp(y(x))+2 = 0:
ics:=y(0)=0.6, D(y)(0)=0:
Sol1:=dsolve({ode,ics}, numeric):
A:=plots:-odeplot(Sol1, [x,y(x)], x=0..5, color=red):
Sol2:=dsolve({ode,ics}, y(x), 'series');
B:=plot(eval(y(x), Sol2), x=0..5, color=blue):
plots:-display(A, B, view=-5..2);

Edit.
 

If I understand correctly  op(0, e)  returns not the type but a procedure (a command, a function) by which this expression  e  is constructed. Sometimes this coincides with a type, but not always. See:

op(0, n!);
op(0, sin(x));
                                                 factorial
                                                    sin

whattype(n!);
whattype(sin(x));
                                                 function
                                                 function

 Integer  is the command in Maple. See help on this. I do not know why this command is needed. I never used it.                                  

Here are all the options to plot a rectangle:

restart;
A:=plots:-display(plottools:-rectangle([0.2, 0.8], [0.8, 0.2], style=line, color=green, thickness=3)):  # Only the border
B:=plots:-display(plottools:-rectangle([0.2, 0.8], [0.8, 0.2], color=green, linestyle=2)):  # Border in the desired style
C:=plots:-inequal({x>0.2, x<0.8, y>0.2, y<0.8}, x=0..1, y=0..1, color=green, nolines):  # Without any border
E:=plots:-display(plottools:-rectangle([0.2, 0.8], [0.8, 0.2]), style=line,  color=red, thickness=3):  # Border in the desired color
plots:-display(< <A | B >, <C | plots:-display(E,C)> >, view=[0..1,0..1], scaling=constrained, axes=none, size=[250,250]);  # "Array" of all the variants

acer's and vv's ways (S  and  s  procedures) don't work explicitly for  n>50  (don't return an explicit integer):

  S(100);
  s(100);
                                          100*exp(1)*GAMMA(100, 1)
                                          100*exp(1)*GAMMA(100, 1)

Here is another way:

r := rsolve({u(n)=n*(u(n-1)+1), u(0)=0}, u):
SS:=s->simplify(value(eval(r, n=s))):
SS(100);

253686955560127297415270748212280220445147578566298142232775185987449253908386446518940485425152049793267407732328003493609513499849694176709764490323163992000

Of course, it's easy to combine with the previous results, so that the explicit result is returned only at the user's request (for those who are annoyed by long numbers).