The domain of definition of your function is a semicircle. Use the direct plotting in this domain:

plot3d(4*x*y^2-x^2, x=0..1, y=-sqrt(1-x^2)..sqrt(1-x^2), scaling=constrained, numpoints=10000, axes=normal);

Another option is the using of the polar coordinates:

plot3d(eval([x,y,4*x*y^2-x^2],[x=r*cos(t),y=r*sin(t)]), r=0..1, t=-Pi/2..Pi/2, scaling=constrained, axes=normal);

VerticesOfCuboid:=proc(V1::list, V2::list)
[V1, seq(subsop(i=V2[i], V1), i=1..3), V2, seq(subsop(j=V1[j], V2), j=1..3)];
end proc: 

Example of use:

VerticesOfCuboid([1,1,2], [3,4,5]);

               [[1, 1, 2], [3, 1, 2], [1, 4, 2], [1, 1, 5], [3, 4, 5], [1, 4, 5], [3, 1, 5], [3, 4, 2]]

Addition. In geom3d you can specify a cuboid by  parallelepiped  command and find its parameters by other commands (detail, vertices, faces  and so on).

There are several ways to do this. Here are two ones:

restart;
f := unapply(x*(8*x^2+5*x+cos(y)), x, y);
p := solve({diff(f(x, y), x) = 0, diff(f(x, y), y) = 0}, {x, y});
eval([x,y], p[1]);  # The first way
assign(p[1]);  # The second way 
x, y; 

I understood  this  to mean that  P  is an operator in the functional space and therefore we can do so:

P:=u->(i,j)->u(i+1,j);

Example of use:

u:=(i,j)->i^2+j^2;
v:=P(u);  # v is a new function
v(2,3);
                                    
 

For  Q  everything is similar.

Example:

y:=t-><sin(t), cos(t), exp(2*t)>;
series~(y(t), t, 3);

For the example, I have specified values for  u, L, c  :

restart;
z:=(-1/(4*(-u+L)))*(4*p^2*b^2*d-3*g*c*p^2);
u:=2:  L:=1:  c:=0.6: 
constr:={ 1-(2*d)/(p*(sqrt(g)+1))<=b, b<=1, 0<=d, d<=L, c<=p, p<=u, 1<=g, g<=4};  
Optimization:-Maximize(z, constr);

Edit.

You can use  `A=B=C`  instead. Output will be without any quotes.

If you need to use this as a logical condition (for example inside a selection statement  if ... fi), then write  A=B and B=C

Try (only for purely visual purposes)

exp(`2 t`) ;

This is a workaround for Classic Interface. In Standard one there are no any parentheses around the exponent.

Edit.

Use  InertForm  package.

Example:

A:=<1,2; 3,4>:
B:=<3,4; 1,2>:
with(InertForm):
L := Parse("A*B"):
Display(L) = A.B;

A+`&lambda;I`=<1,2; 3,4>;  # Or

A+`&lambda; I`=<1,2; 3,4>;                     

To solve this problem,  combinat:-composition  command is useful. 

restart;
mondeg:=proc(var::list(name), d::nonnegint)
local n, L;
uses combinat, ListTools;
n:=nops(var);
L:=[seq(map(t->t-~1, Reverse(convert(composition(k+n, n),list)))[], k=0..d)];
map(t->`*`((var^~(t))[]), L);
end proc:

Example of use:

mondeg([x,y,z], 3);
         [1, x, y, z, x^2, x*y, x*z, y^2, y*z, z^2, x^3, x^2*y, x^2*z, x*y^2, x*y*z, x*z^2, y^3, y^2*z, y*z^2, z^3]

P:=6*x^2+x-2:
factor(P);
map(t->lcoeff(t)*``(t/lcoeff(t)), factor(P));
                                                             (3*x+2)*(2*x-1)
                                                            6*(x+2/3)*(x-1/2)

Does this match what you wanted?

For each n, I took 5 frames to lengthen the whole process:

restart;
with(plots):
f := x -> x;
L := 1: 
fs:=n->add(2*(-1)^(1+i)*sin(Pi*i*x)/Pi/i, i=1..n):
F:=n->plot([fs(n),f(x)], x=-L..L, color=[blue,red]):
display(seq(F(n)$5, n=1..25), insequence);

      


An alternative way to create an animation is to use  plots:-animate  command. Here is the code for your example:

restart;
with(plots):
f := x -> x;
L := 1: 
fs:=n->add(2*(-1)^(1+i)*sin(Pi*i*x)/Pi/i, i=1..n):
F:=n->plot([fs(n),f(x)], x=-L..L, color=[blue,red]):
animate(display,['F'(floor(n))], n=1..25, frames=100);

This method is preferable if the animation parameter changes in some continuous interval  a .. b .

Edit.

1. In the symbolic solution of a cubic equation, there is usually the imaginary unit  I, even if all the roots are real numbers (the so-called irreducible case). For example:

solve(x^3-3*x+1=0);
fsolve(x^3-3*x+1=0);

2. In general case, it is unlikely that something will be simplified. In special cases in addition simplify command you can also try  radnormal command.

3. In fact, this is the real number:
evalf(0.1159017153e-1*(820.0133-(3168.172220*I)*sqrt(3))^(1/3)+3.632452480/(820.0133-(3168.172220*I)*sqrt(3))^(1/3)+.6560037088+(1/2*I)*sqrt(3)*(-0.2318034307e-1*(820.0133-(3168.172220*I)*sqrt(3))^(1/3)+7.264904960/(820.0133-(3168.172220*I)*sqrt(3))^(1/3)));

                       0.6965691596-5.660254040*10^(-11)*I

The imaginary unit appears as a result of rounding errors.

In any case, you somehow have to refer to these rows. This is convenient if you use end-to-end indexing.

For your example:

row[1]:=[1,2,3];  
row[2]:=[4,5,6,7];
row[3]:=['red','blue',7];
Matrix(convert(row, list));

Note that Maple automatically aligns the lengths of the rows by adding zeros.

If you want to keep empty spaces in case of different lengths of lists, you can do this:
row[1]:=[1,2,3]:  
row[2]:=[4,5,6,7]:
row[3]:=['red','blue',7]:
m:=max(seq(nops(row[k]), k=1..3));
Matrix([seq([row[i][], ``$(m-nops(row[i]))], i=1..3)]);

Edit.