Here is an example of animation of matrixplot . It was exported as  gif - file:

A:=LinearAlgebra:-RandomMatrix(7, generator = 0 .. 9);
P:=plots:-matrixplot(A, heights = histogram, axes = normal):
f:=t->plottools:-rotate(P, t, [[0,0,0],[0,0,1]]):
plots:-animate(f,[t], t=0..2*Pi, frames=90, scaling=constrained); 

Addition. Another example - the rotation with the axes:

A:=LinearAlgebra:-RandomMatrix(7, generator = 0 .. 9):
plots:-matrixplot(A, heights = histogram, axes = normal):
f:=t->plots:-matrixplot(A, heights = histogram, axes = normal, orientation=[-t,60]):
plots:-animate(f,[t], t=0..360, frames=90, scaling=constrained);

ThirdDer := proc(f::procedure, x0, h)
(f(x0+2*h)-2*f(x0 +h)+2*f(x0 -h)-f(x0 -2*h))/2/h^3;
end proc: 

Example of use:

ThirdDer(sin, Pi, 0.01);

                                               0.9999800000

See these results in Google search

https://www.google.ru/?gws_rd=ssl#newwindow=1&q=SADE+package+for+Maple&*

extrema command does not find the local maxima or minima of a function, but only finds the values of the function at critical points (and these points themselves with  's'  option), i.e. the points at which all partial derivatives are equal to 0. Among these points there can be points of a maximum or a minimum or saddle points. Therefore, all points found by  extrema  command must be further investigated using the values of the second partial derivatives at these points. In Maple, this is automated by  Student:-MultivariateCalculus:-SecondDerivativeTest  command. See the solution of the first example:

restart;
extrema(x^3+y^3-3*x*y, { }, {x,y}, 's');  # The values of the function at the points  [1,1]  and  [0,0]
s;  # The critical points themselves 
Student:-MultivariateCalculus:-SecondDerivativeTest(x^3+y^3-3*x*y, [x,y]=[0,0]);
Student:-MultivariateCalculus:-SecondDerivativeTest(x^3+y^3-3*x*y, [x,y]=[1,1]);
plot3d(x^3+y^3-3*x*y, x=-1..2.7, y=-1..2.7, color=khaki, view=-2..1.7, axes=normal);

Thus we see that  only at the point  [1,1]  the local minimum is reached (the minimum is -1), and the point  [0,0]  is a saddle point. The plot confirms these conclusions.

This method  (with  Student:-MultivariateCalculus:-SecondDerivativeTest  command) is not suitable for solving the second example, because here there is the constraint equation. Fortunately, the constraint equation  x²+y²=5  defines a compact set in the plane (this is a circle) on which the continuous function  z=x+2*y  takes its smallest and largest values. But since we get only 2 critical points, one of these points will be the maximum, and the other will be the minimum:

Expr:=x+2*y:  constr:=x^2+y^2=5:
extrema(Expr, constr, {x,y}, 's');
s;

Thus we see that  at the point  [-1, -2]  the local minimum is reached (the minimum is  - 5),  and at the point  [1, 2]  the local maximum is reached (the maximum is  5) .

plot(sin(x), x=-Pi..Pi, style=point, color="Red", legend="sin(x)", symbol=diamond, numpoints=20, adaptive=false, scaling=constrained);

Here is a numerical calculation of this double sum:

Download CoilCalculus_new.mw

You just did not simplify the result:

restart;
util:=x1^(1/3)*x2^(1/2):
constraint:=m-P1*x1-P2*x2:
maxutil:=util+lambda*constraint;
CP01:=diff(maxutil,x1);
CP02:=diff(maxutil,x2);
CP03:=diff(maxutil,lambda);
solset:=solve({CP01=0,CP02=0,CP03=0}, {lambda,x1,x2}, explicit);
simplify(solset[1]) assuming m>0, P1>0, P2>0;
 

restart;
M:=Matrix([[a,b,c], [d,e,f], [g,h,k]]);
M[1];   # First row
M[..,1];   # First column
M[2..3];   # Second and third rows
M[..,1..2];   # First and second columns
M[..,[1,3]];   # First and third columns
M[1..2, 2..3];   # First and second rows, second and third columns

Addition: 3 more examples:

M^%T;   # Transpose a matrix
M^(-1);   # Inverse matrix
M.M^%T;   # Multiplication of matrices  (using the usual dot)

f:=unapply(simplify(sum(i, i=1..n)), n);
f(n);
f(10);

Unfortunately there are no solutions to this system that completely satisfy your conditions. First we find all the solutions of this system (real and imaginary). Maple easily solves polynomial systems, so first we squared the last 2 equations. The resulting system will be equivalent to the original system under conditions  x>=0, x+y>=0 .

restart;
eq1:=4.1*10^9-x*Tl^4-4.527587408*10^7/Pi=0:
eq2:=4.1*10^9-(0.6013906720-x-y)*Th^4+4.527587408*10^7/Pi=0:
eq3:=4.1*10^9-Ti^4*y=0:
eq4:=1-1.487209994*10^(-20)*Tl^8/(1+sin(-0.41)^2)=x^2:
eq5:=1-1.487209994*10^(-20)*Th^8/(1+sin(-0.41)^2)=(x+y)^2:
_MaxSols:=infinity:
sol:={solve({eq1,eq2,eq3,eq4,eq5}, {x,y,Tl,Th,Ti})}:
sol1:=simplify(sol, zero):
nops(sol1);

                                                              1024

Maple found 1024 solutions, among which there are many imaginary solutions. Next we first try to take all the real solutions and see that there are no such solutions where all the roots are real. But there are solutions where all the roots are real except for Ti :

S1:=select(t->eval(`and`(Im(x)=0, Im(y)=0, Im(Tl)=0, Im(Th)=0, Im(Ti)=0), t), sol1);
S2:=select(t->eval(`and`(Im(x)=0, Im(y)=0, Im(Tl)=0, Im(Th)=0), t), sol1):
nops(S2);
S3:=select(t->eval(`and`(x>=0, x+y>=0, Tl>Th), t), S2);
nops(S3);

This can be done in a mass of ways. If you want the cone to have the specified proportions, you first need to do some calculations to find out the necessary parameters (vertex coordinates, base radius and height).

The first (easy) way:

A:=plottools:-cone([0,0,-1], 1, 1, style=surface, color="Green"):
B:=plottools:-sphere(color=pink, transparency=0.5):
C:=plottools:-reflect(A, [[0,0,0],[1,0,0],[0,1,0]]):
plots:-display(C, B, axes=normal, view=[-1.4..1.4,-1.4..1.4, -1.2..1.4], scaling=constrained);

Addition. The second way.

 Let us  inscribe a cone , whose height is twice the radius of the base, in the unit sphere. First we need to find the radius of the base of this cone:

R:=solve(2*x=sqrt(1-x^2)+1);  # Radius of the base of the cone
A:=plot3d([r*cos(t), r*sin(t), -sqrt(1-R^2)], r=0..R, t=0..2*Pi, style=surface, color="Green"):  # The base of the cone
B:=plot3d([r*cos(t), r*sin(t), -2*r+1], r=0..R, t=0..2*Pi, style=surface, color="Green"): # Lateral surface of the cone
C:=plottools:-sphere(color=pink, transparency=0.7):
plots:-display(A, B, C, axes=normal, view=[-1.2..1.2,-1.2..1.2, -1.2..1.4], scaling=constrained);

Here the cone surfaces are given in terms of their parametric equations.

I supposed that  A, B, alpha, y  are reals,  i=I  is imaginary unit.

restart;
w := A*exp(-alfa*(1+I)*y)+B*exp((1+I)*y);            
u := evalc(w);
subs(I=0, u);

eq:=diff(y(x),x)=sqrt(1+a*x+2*y(x)):
eq1:=eval(eq, a=0);  # Substitution a=0 into the equation eq
sol:=dsolve({eq1, y(0)=1}, y(x));  # Solution of the equation eq1
eval(eq1, sol);  # Substitution of the solution into the equation eq1
is(%) assuming x>0;  # Check of the solution
eval(sol, x=1);  # Exact value y(1)
evalf[5](%);  # Approximate value y(1)

F:=unapply(int(1/(0.1060444429e-1-0.2120888857e-1*X+0.1033933318e-1*X^2), X), X);
fsolve(F(x)-F(0)=5000, x);
simplify(fnormal(%), zero);

See  LinearAlgebra:-GenerateMatrix  command in Help.