Another way is to use  RealDomain  package:
 

with(RealDomain):
plot(x^(5/3)-5*x^(2/3), x=-3..7, color=red);

In the first example, we use the standard parametrization of the sphere and a rotation. In the second example (the simpler one) the standard parametrization only.

The first example (the upper and lower parts after cutting):

P1:=plot3d([4*cos(phi)*sin(theta),4*sin(phi)*sin(theta),4*cos(theta)], phi=0..2*Pi, theta=0..arccos(-1/2/sqrt(3)), style=surface, color=khaki):
P2:=plot3d([4*cos(phi)*sin(theta),4*sin(phi)*sin(theta),4*cos(theta)], phi=0..2*Pi, theta=arccos(-1/2/sqrt(3))..Pi, style=surface, color=khaki):
P:=plot3d([r*cos(phi),r*sin(phi),-2/sqrt(3)], phi=0..2*Pi, r=0..sqrt(44/3), style=surface, color=khaki):
S1:=plottools[rotate](P1, arccos(1/sqrt(3)), [[0,0,0],[1,-1,0]]):
S2:=plottools[rotate](P2, arccos(1/sqrt(3)), [[0,0,0],[1,-1,0]]):
S:=plottools[rotate](P, arccos(1/sqrt(3)), [[0,0,0],[1,-1,0]]):
plots[display](<plots[display](S1,S) | plots[display](S2,S)>, axes=normal, view=[-5.9..5.9,-5.9..5.9,-4.9..4.9], orientation=[55, 75, 0], scaling=constrained);

The second example (The left and right parts after cutting):

P1:=plot3d([4*cos(theta),4*cos(phi)*sin(theta),4*sin(phi)*sin(theta)], phi=0..2*Pi, theta=0..arccos(5/16), style=surface, color=khaki):
P2:=plot3d([4*cos(theta),4*cos(phi)*sin(theta),4*sin(phi)*sin(theta)], phi=0..2*Pi, theta=arccos(5/16)..Pi, style=surface, color=khaki):
P:=plot3d([5*cos(theta)-2,5*cos(phi)*sin(theta),5*sin(phi)*sin(theta)], phi=0..2*Pi, theta=0..arccos(13/20), style=surface, color=khaki):
plots[display](<plots[display](P1,P) | plots[display](P2,P)>, axes=normal, view=[-5.9..5.9,-5.9..5.9,-4.9..4.9], orientation=[55, 75, 0], scaling=constrained);

All 3 methods are incorrect. Here is the corrected first method:

T := plottools[transform]((x, y, z) ->`if`(z <= 2-y, [x, y, z], [x, 1+y/2-z/2, 1-y/2+z/2])):
plots[display](T(plottools[sphere]([0, 0, 0], 4, style=surface, color=khaki, lightmodel=light4, numpoints = 100000)), scaling = constrained);

Clarification. [x, 1+y/2-z/2, 1-y/2+z/2]  is the projection of the point  [x, y, z]  onto the plane  z = 2 - y

Edit.

func := n -> map(t -> t/content(t), {op(expand((a+b+c)^n))}):

Example of use:
func(3);

                     {a^3, b^3, c^3, a*b^2, a*c^2, a^2*b, a^2*c, b*c^2, b^2*c, a*b*c}

If you have a planar region, use the double integral. See this article

If the boundary of a plane region is specified as a piecewise smooth curve, then to calculate double integrals is advisable to use Green's formula. This technique I've used here to calculate the area of a plane figure.

Try  Student[MultivariateCalculus][LagrangeMultipliers]  command with  output=detailed  option.

For numurical solution use  Optimization[Maximize]  or   Optimization[Minimize]  commands.

To use  plot3d  command, you must from each plane equation  Plane 1 .. Plane 3  express  z  as a function of  x  and  y :

restart;
P1:= x + 3*y - 5*z = 0:
P2:= x + 4*y - 8*z = 0:
P3:= - 2*x - 7*y + 13*z = 0:
S:=map(solve, {P1,P2,P3}, z);
plot3d(S, x = -8 .. 8, y = -10 .. 10,  style=surface, color = [blue,red,green], axes=normal); 

Use formulae of change from spherical coordinates  [r,theta,phi]  to cartesians coordinates  [x,y,z] :                 

x = r*cos(phi)* sin(theta)
y = r*cos(phi)* sin(theta)
z = r*cos(theta)

Example of use:

S:=[2, Pi/4, Pi/3]: # Spherical coordinates
L:=[r, theta, phi]=~S;
CartCoord:= [r*cos(phi)*sin(theta), r*sin(phi)*sin(theta), r*cos(theta)]:
eval(CartCoord, L);  # Cartesian coordinates
 

Use nested  seq  command for Q5 and Q6 questions.

For example for Q6:

F[6] = {seq(seq(p/q, p=1..q), q=1..6)};

If I understand your question then

A := [0., 12.0, 37.0, 62.0, 87.0, 112.0, 137.0, 162.0, 187.0, 212.0, 237.0, 238.0]:
B := [8.0, 8.0, 7.0, 6.0, 5.0, 4.0, 3.0, 2.0, 2.0, 1.0, 1.0, 0.]:
n:=nops(A):
B:=unapply(piecewise(t>=A[1] and t<=A[2],B[1], seq(op([t>A[i] and t<=A[i+1],B[i+1]]), i=2..n-1), t>238,0), t);
B(30), B(60);
plot(B(t), t=0..250, discont);

The simple procedure  LMsystem  solves your problem. The required formal parameters of the procedure: expr  is the objective function specified as an expression, constr  is the list of constraint functions, assumed equal to 0, vars  is list of names of independent variables. The optional parameter  opt  (to be traditional (by default) or set or list) specifies the form of output for the system of equations.

LMsystem:=proc(expr, constr::list, vars::list, opt:=traditional)
local n, m, F, lambda, L;
n:=nops(vars);
m:=nops(constr);
F:=expr+add(lambda[i]*constr[i], i=1..m);
L:=[seq(diff(F,vars[i])=0, i=1..n), seq(constr[j]=0, j=1..m)];
if opt=traditional then %piecewise(seq(op([``, L[i]]), i = 1 .. nops(L))) else
if opt=set then {L[]} else L fi; fi;
end proc:

Examples of use:

LMsystem(2*x+y^2-2*z^2, [2*x+y-5*z-6, -x+3*y-z-3], [x,y,z]);  # Output in a traditional form

LMsystem(2*x+y^2-2*z^2, [2*x+y-5*z-6, -x+3*y-z-3], [x,y,z], set);  # Output as a set

If you use Joe's hint:

f1:=(1/s^2)^N:
simplify(f1)  assuming s>0;

                                               s^(-2*N)

LinearAlgebra[Eigenvalues]  command returns the column vector whose entries are all the eigenvalues of the given matrix.

Example of plotting:

restart;
A:=LinearAlgebra[RandomMatrix](4, datatype = float);
V:=LinearAlgebra[Eigenvalues](A);
z1, z2, z3, z4:=seq(V[i], i=1..4);
c:=plot([Re,Im]~([%]), style=point, symbol=solidcircle, color=red, symbolsize=15):
T:=plots[textplot]([seq([(Re,Im)(z||i),cat('z',i)], i=1..4)], align={right,above}, font=[times,roman,16], color=red): 
plots[display](c, T, scaling=constrained);

Edit.

From wiki: "Two matrices are equivalent if and only if they have the same rank"

So use  LinearAlgebra[Rank]  command.

You can easily put labels on those places where you want, if you use  plots[textplot3d]  command. The function must be specified as a procedure, rather than the expression (see Carl's comments below).

Example:

S:=plot3d((x,y)->x^2+y^2, -2.4..2.4, -2.4..2.4, axes=normal):
Labels:=plots[textplot3d]([[2.4,0,0,"x", align=above],[0,2.4,0,"y", align=above], [0,0,12,"z", align=left]], font=[times,roman,16]):
plots[display](S, Labels, axes=normal);