Probably the paraboloid must be round and 0.25 is the maximum z-coordinate.

restart;
solve({eval(x^2/2+y^2/2, [x=R,y=0])=0.25, R>0});
assign(%):
plot3d([r*cos(phi),r*sin(phi),r^2/2], r=0..R, phi=0..2*Pi, scaling=constrained, axes=normal);

For the complete solution of the problem, you need to find the surface area of this paraboloid.

In  GraphTheory  package there is no  copy  command, but there is  CopyGraph  command.

Addition. Example:

with(GraphTheory):
G := CycleGraph(5);
H := CopyGraph(G);
V:=[1,2,3]; 
HighlightVertex(H, V, red);
 DrawGraph(H);
 DrawGraph(G);
 

solve(20*x^3+10*x^2+4*x+1, explicit):  
20*(x-%[1])*map(simplify, sort(collect((x-%[2])*(x-%[3]), x)));   # Exact factorization
evalf(%);   # Approximate factorization

Let  x=a*t ,  y=b*t+c

We have  cos(x)=cos(y)  and  sin(x)=sin(y)

These equalities are equivalent  to   x - y = 2*Pi*k ,  k  is integer

So  t = (c+2*Pi*k)/(a-b)

Conclusion: thus we see that the set of all solutions is the infinite arithmetic progression with initial term  c/(a-b)  (if  k=0)  and the common difference  2*Pi/(a-b)

Addition. The example of finding all solutions in a given interval  (0..20) for specific parameter values  [a=7,b=2,c=3]:

R:=eval((c+2*Pi*k)/(a-b), [a=7,b=2,c=3]);
K:=[isolve({R>=0,R<=20})];
seq(eval(R,k), k=K);

The first procedure uses  Partition  procedure from here. It gives all possible partitions of an odd number  >=9  into 3 prime summands.

PartitionOfOdd:=proc(n::odd)
local p;
uses NumberTheory;
p:=prevprime(n);
Partition(n, 3, {seq(ithprime(k),k=2..pi(n))});
end proc:

Example of use:

PartitionOfOdd(63);

[ [3, 7, 53], [3, 13, 47], [3, 17, 43], [3, 19, 41], [3, 23, 37], [3, 29, 31], [5, 5, 53], [5, 11, 47], [5, 17, 41], [5, 29, 29], [7, 13, 43], [7, 19, 37], [11, 11, 41], [11, 23, 29], [13, 13, 37], [13, 19, 31], [17, 17, 29], [17, 23, 23] ]

The second procedure gives only 1 partition, but it is much faster than the first one:

PartitionOfOdd1:=proc(n::odd)
local N, i, j, a, b, c;
uses NumberTheory;
N:=pi(n);
for i from 2 to N-1 do
for j from 2 to N-1 do
a:=ithprime(i); b:=ithprime(j); c:=n-a-b;
if isprime(c) then return [a,b,c] fi;
od; od;
end proc: 

Example of use:

PartitionOfOdd1(123456789);

                                                           [3, 29, 123456757]

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);