To plot one or more points, the simplest way will be to use the syntax  plot([A, B, ...], style=point) , where A, B, ...  are points given by lists of their coordinates. To specify the color, size and shape of symbols representing points, use the options  color =... ,  size = ...  and  symbol = ...  . To label points we use  plots:-textplot  command. Note that the function  y=x^2  is defined as a procedure  f:=x->x^2 . This is convenient when calculating the values of both the function itself and its derivative at specific points.

Below, as an example, we solve the problem:

Given the curve  y=x^2  and 2 points  A(-1,1) ,  B(2,4)  .
1. Check that these points lie on this curve.
2. Find the point on the curve  y=x^2 and the equation of the tangent line to this curve (at this point)  such that the tangent is parallel to the chord  AB .
3. Draw all the objects in one plot.

Solution:

restart;
f:=x->x^2:  A:=[-1,1]:  B:=[2,4]:

# Calculations
#
x1, x2, y1, y2 := A[1], B[1], A[2], B[2];
is(f(x1)=y1), is(f(x2)=y2); # Checking that points A and B lie on the curve y=x^2
Eq := D(f)(x0)=(f(x2)-f(x1))/(x2-x1);  # Mean value theorem
x0:=solve(Eq, x0);
y0:=f(x0);
k:=D(f)(x0); # The slope of the tangent
Tangent:=y=y0+k*(x-x0); # The equation of the tangent 
C:=[x0,y0];  # The point of tangency

# Plotting
#
Lines:=plot([f(x), [A,B], rhs(Tangent)], x=x1-1..x2+0.5, color=[blue,green,red]): # The plot of all the lines
Points:=plot([A,B,C], style=point, color=blue, symbol=solidcircle, symbolsize=12): # The plot of all the points
Labels:=plots:-textplot([[A[],"A"],[B[],"B"],[C[],"C"],[-1,3.5,y=f(x)],[-1,-1.5,Tangent]], font=[times,16], align={left,above}): # All the labels
plots:-display(Lines,Points,Labels, scaling=constrained, view=[x1-1..x2+0.5,-2..5], size=[500,500]);

When you apply thу  simplify  command to some equality, it is separately applied to the left and right sides, i.e.  simplify(A=B)  is equivalent to  simplify(A)=simplify(B)  (that's the design). 
To check some equality (an identity) use the  is  command instead of simplify:

restart;
simplify((a^2-b^2)/(a-b)=a+b);
is((a^2-b^2)/(a-b)=a+b);

If you really want to simplify your expression  f  to the expression  g , there is a workaround. We make simple substitutions, after which the numerator and denominator will be polynomials, then simplify. I wonder if the Mathematica can directly simplify  f  to  g  without such tricks.

restart;
f:=(9*(x^(-2/3*a))^2*exp(6/a*(x^(-2/3*a))^(1/2))^2*_C0^2-6*(x^(-2/3*a))^(3/2)*exp(
6/a*(x^(-2/3*a))^(1/2))^2*_C0^2*a+x^(-2/3*a)*exp(6/a*(x^(-2/3*a))^(1/2))^2*_C0^
2*a^2+18*(x^(-2/3*a))^2*exp(6/a*(x^(-2/3*a))^(1/2))*_C0-2*x^(-2/3*a)*exp(6/a*(x
^(-2/3*a))^(1/2))*_C0*a^2+6*(x^(-2/3*a))^(3/2)*a+x^(-2/3*a)*a^2+9*(x^(-2/3*a))^
2)/(3*_C0*exp(6/a*(x^(-2/3*a))^(1/2))*(x^(-2/3*a))^(1/2)-exp(6/a*(x^(-2/3*a))^(
1/2))*_C0*a+3*(x^(-2/3*a))^(1/2)+a)^2:
g:=x^(-2/3*a):
Subs:=[x^(-2*a*(1/3))=A^2,exp(6*sqrt(x^(-2*a*(1/3)))/a)=B];
f1:=subs(Subs, f);
f2:=simplify(f1) assuming A>0;
subs(A=sqrt(x^(-2*a*(1/3))), f2);
is(%=g);

restart;
A:=Matrix(3, [$1..9]);
L:=[indices(A, pairs)];
L1:=(rhs=lhs)~(L);
T:=table(L1);
T[7];

Eq:=a^(5/2)*s^2*(D(f))(eta)^2*sqrt(nu)/(R*sqrt(a)*sqrt(nu)+eta*nu) = a^(5/2)*s^2*(D(P(s, eta*sqrt(nu)/sqrt(a))))(eta)/sqrt(nu);
Term:=select(has,rhs(Eq),P);
solve(Eq, Term);

It's better to implement  this as a procedure so that  Digits  is set automatically:

restart;
Evalf:=proc(x::realcons,n::posint)
local m;
m:=ilog10(evalf(x));
Digits:=m+n+3;
printf(cat("%.",n,f),x);
end proc:

Examples of use:

Evalf(123456789/123,6);
``;
Evalf(123456789/123,15);
``;
Evalf(Pi,3);

If I understood everything correctly, then your surfaces can be plotted according to the following scheme:
1. We define the parametric equations of the curve:  x = R*cos(phi),  y = a - R*(1 - sin(phi)) (phi is the parameter) in the plane  xOy , setting the corresponding constants  a  and  R .
2. Both surfaces are obtained by rotating these curves around the axis  Ox .

In the picture in 3D below, the orientation of the axes  x, y, z  roughly matches your original drawing ( x  is down, y  is back, z  is right):

restart;
a:=3: R:=8:  # For the curve A 
a1:=2: R1:=-8:  # For the curve B
x:=R*cos(phi); y:=a-R*(1-sin(phi)); x1:=R1*cos(phi); y1:=a1-R1*(1-sin(phi));
A:=plot([x, y, phi=Pi/3..2*Pi/3], title="The plot A"):
B:=plot([x1, y1, phi=Pi/3..2*Pi/3], title="The plot B"):
plots:-display(< A | B >, scaling=constrained, view=[-4..4,-4..4], labels=["x","y"], labelfont=[times,bold,16], titlefont=[times,16]);

A1:=plot3d([x,y*cos(alpha),y*sin(alpha)], phi=Pi/3..2*Pi/3, alpha=0..2*Pi, axes=normal, title="The plot A1"):
B1:=plot3d([x1,y1*cos(alpha),y1*sin(alpha)], phi=Pi/3..2*Pi/3, alpha=0..2*Pi, axes=normal, title="The plot B1"):
plots:-display(< A1 | B1 >, scaling=constrained, view=[-4.5..4.5,-4.5..4.5,-4.5..6.5], labels=["x","y","z"], labelfont=[times,bold,16], titlefont=[times,16],  orientation=[-10,-3,70]);

As an addition - animation of  A1 :

plots:-animate(plot3d,[[x,y*cos(alpha),y*sin(alpha)], phi=Pi/3..2*Pi/3, alpha=0..'a', axes=normal], 'a'=0..-2*Pi, frames=90, scaling=constrained, view=[-4.5..4.5,-4.5..4.5,-4.5..4.5], labels=["x","y","z"], labelfont=[times,bold,16], orientation=[-10,-3,70]);

You can use the  subs  command if you notice that  Zs / Z_AB = K_S  is equivalent to  Zs=Z_AB*K_S :

restart;
eq_5 := x = -(K_i*Zs - Z_AB - Zs)/(K_i*Z_AB);
eq_5_m5 := x = expand(rhs(eq_5));
eq_5_m6 := x = subs(Zs=Z_AB*K_S, rhs(eq_5_m5 ));

I usually use the following way:

restart;
u:=<-2+m,3+m>;

u:=eval(u, m=5);
# or
u1:=eval(u, m=5);

Such equations are easily reducible to ordinary differential equations with derivatives. In the example below, we solve the equation  p*dx-q*dy=1/2*d(y^2-x^2)  by taking  p=3, q=2  and plotting several integral curves and 2 asymptotes:

restart;
d:=z->diff(z,x)*dx:
p:=3: q:=2:
p*d(y(x))-q*dx=1/2*d(y(x)^2-x^2);
dsolve(%);
plot([seq(eval(y(x),%[1]),_C1=-8..3),seq(eval(y(x),%[2]),_C1=-8..3),-(x-2)+3,x-2+3], x=-4..8, color=[green$24,red$2], scaling=constrained);

restart;
S := proc(n) 
local i, s;
uses NumberTheory;
s:=0; 
for i to n do 
if i in Divisors(n) and type(n/i, odd) then s:=s+i end if; 
end do;
s;
end proc:

# Example
S(12);

                                                  16

The code below can be made more automatic, but it's a rather tedious job. Let it be for now:

restart;
expr1:=sqrt(x^2*y - 4)*sqrt(x^2*y);
expr2:=sqrt((x^2*y - 4)*(x^2*y));
SolveTools:-SemiAlgebraic({op([1,1],expr1)>=0,op([2,1],expr1)>=0, op([1,1],expr2)>=0 , expr1^2=expr2^2}, {x,y});

restart;

plot3d(y^2, x=0..1, y=-1..1, style=surface, color=grey, filled=true, axes=normal, view=[-0.5..1.5,-1.5..1.5,0..1.5]);

Volume=int(y^2, [y=-1..1,x=0..1]);

plots:-implicitplot(x=y^2, x=0..9, y=-3..3);

The programm uses the first  2*K terms of your list only (you have  K=1 ). So all works properly:

restart;
interface(version);

genfunc:-rgf_findrecur(1, [1, 2], t, n);
t:=unapply(rsolve({%, t(1)=1},t(n)), n);
seq(t(n), n=1..4);

The best quality will be if all surfaces are parameterized:

restart;
x:=r*cos(phi); y:=r*sin(phi); z:=3-x; H:=eval(z,r=2);
S1:=plot3d([x,y,0], r=0..2, phi=0..2*Pi, style=surface, color=blue): # Bottom
S2:=plot3d([eval([x,y],r=2)[],h], phi=0..2*Pi, h=0..H, style=surface, color=green): # Side surface
S3:=plot3d([x,y,z], r=0..2, phi=0..2*Pi, style=surface, color=red): # Top surface
plots:-display(S1,S2,S3, axes=normal, view=[-3..3,-3..3,0..5], scaling=constrained, lightmodel=light4, labels=["x","y","z"], labelfont=[times,bold,14]);

We can also show the cutting plane using the  transparency  option:

restart;
x:=r*cos(phi); y:=r*sin(phi); z:=3-x; H:=eval(z,r=2);
S1:=plot3d([x,y,0], r=0..2, phi=0..2*Pi, style=surface, color=blue):
S2:=plot3d([eval([x,y],r=2)[],h], phi=0..2*Pi, h=0..H, style=surface, color=green):
S3:=plot3d([x,y,z], r=0..2, phi=0..2*Pi,  style=surface, color=red):
S4:=plot3d(3-u, u=-3..3, v=-3..3, style=surface, color=yellow, transparency=0.5):
plots:-display(S1,S2,S3,S4, axes=normal, view=[-3..3,-3..3,0..5], scaling=constrained, lightmodel=light4, labels=["x","y","z"], labelfont=[times,bold,14]);

                          

Edit.