If I understand your problem  you can use  seq  command.

Example of use:

seq(LinearAlgebra[RandomMatrix](3), i=1..3);

A variation on the theme vv's solution:

restart; 
S := {-x+2*y, -x-y+3*z, 3*x-2*y-z, x+y+z = 1}: 
A, B, C:=map2(eval, [x,y,z], solve~([seq(S minus {S[i]}, i = 1 .. 3)]))[ ]; 
OO := [0, 0, 0]: 
plots:-polygonplot3d([[OO, A, B], [OO, A, C], [OO, B, C]], axes = normal, labels = [x, y, z], labelfont = [times, roman, 16]);

Solid_angle.mw

Edit.

Many integrals (even without parameters ) can not be calculated exactly (symbolically). In this case, the numerical calculation can be used with the specific values of the parameters:

f := GAMMA(phi)*y^(mu*phi-1)*(1-y)^((1-mu)*phi-1)/(GAMMA(mu*phi)*GAMMA((1-mu)*phi)); 
int(eval(ln(1-y)*f, {mu = 1/3, phi = 2}), y = 0 .. 1); 
evalf(%);

v1:=r->rhs(TOV[5])(r)-(rhs(TOV[5])(608529.7)+0.36555);  # The new function
v1(100), v1(1000), v1(100000);  # Examples of use
plot(v1, 0..600000);  # The plot

PDE := {diff(phi(x, t), t) = (diff(phi(x, t), x, x))/S__c-K__r*phi(x, t)+S__r*(diff(theta(x, t), x, x)), diff(u(x, t), t) = diff(u(x, t), x, x)-M^2*(u(x, t)-m*w(x, t))/(m^2+1)-u(x, t)/`ϰ`-2*Omega^2*w(x, t)+Gr*theta(x, t)+Gm*phi(x, t), diff(w(x, t), t) = diff(w(x, t), x, x)+M^2*(m*u(x, t)-w(x, t))/(m^2+1)-w(x, t)/`ϰ`+2*Omega^2*u(x, t), diff(theta(x, t), t) = lambda*(diff(theta(x, t), x, x))/P__r}:
IBC:={phi(0, t) = 1, phi(9, t) = 0, phi(x, 0) = 0, u(0, t) = t, u(9, t) = 0, u(x, 0) = 0, w(0, t) = 0, w(9, t) = 0, w(x, 0) = 0, theta(0, t) = 1, theta(9, t) = 0, theta(x, 0) = 0}:
P__r:=0.71: lambda:=1.0: S__c:=0.22: K__r:=1.0: S__r:=0.5: m:=0.5: M:=sqrt(10.0): `ϰ`:=0.5: Omega:=sqrt(5.0): Gr:=6.0: Gm:=5.0:
sol:=pdsolve(PDE , IBC, numeric, spacestep=0.01);
p1 := sol:-plot(t = 0.3, color = red):
p2 := sol:-plot(t = 0.5, color = gold):
p3 := sol:-plot(t = 0.7, color = yellow):
p4 := sol:-plot(t = 1., color = green):
plots[display]({p1, p2, p3, p4}); 

Edit: spacestep  option has been reduced.

Addition.  The plots of the function phi(x,t) (by default the first function of the system) have been plotted above for the specific values of  t .  If you need the plots of the other functions of the system (for example  u(x,t) ) , you can do it like this:

q1, q2, q3, q4:=seq(eval(u(x,t), sol:-value(t = t0, output = listprocedure)), t0=[0.3, 0.5, 0.7, 1]):
plot([q1,q2,q3,q4], 0..3, color=[red,gold,yellow,green], scaling=constrained, labels=[x, u(x,t)]);

Should be  Md  instead of  Mb :

solve({Fh, Fz, Md, Ms1, Ms2}, {Az, Bz, Dx, Dz, Ma});

      {Az = 44.-1.*Bz-1.*Dz, Bz = Bz, Dx = 1.500000000, Dz = Dz, Ma = 414.-6.*Bz-12.*Dz}

We see that the system has an infinite number of solutions that depend on two parameters  Bz  and  Dz

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