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.

f:=1/(z^2+I)^2;
convert(f, parfrac, sqrt(2));

We obtain this number  sqrt(2)  (a field expansion) by solving the equation

solve(z^2+I, z);

Some expressions are written ambiguously. It is not clear whether the expressions  2*x+x^2  and  sqrt(4*x^3+x^4)  must be in the numerator or in the denominator. I assumed that in the numerator.

restart;
P1:=plots:-shadebetween(0, 1/2*(2*x+x^2)+1/2*sqrt(4*x^3+x^4), x=0..1/13*(5-2*sqrt(3)), color="LightBlue"):
P2:=plots:-shadebetween(1/2*(2*x+x^2)+1/2*sqrt(4*x^3+x^4), 2*x,x=0..1/13*(5-2*sqrt(3)), color="Pink"):
T:=plots:-textplot([[0.1,0.05,c1],[0.1,0.17,c2]], font=[times,bold,16]):
plots:-display(P1,P2,T, view=[0..0.15,0..0.25], labels=[x,z]);

Edit.

restart;
ode:=x*diff(y(x),x)*y(x) = (y(x)^2-9)^(1/2):
ic:=y(exp(4)) = 5:
sol:=dsolve([ode,ic],y(x));
res:=odetest(sol,ode);
solve(res) assuming x>0;

Unfortunately, assuming real  often does not work and a more specific assumption is needed.

Edit.  Below is another solution, in which we first automatically find the domain of the expression  res , and then use this in  assuming :

restart;
ode:=x*diff(y(x),x)*y(x) = (y(x)^2-9)^(1/2):
ic:=y(exp(4)) = 5:
sol:=dsolve([ode,ic],y(x));
res:=odetest(sol,ode);
F:=indets(res, function);
R:=[solve(`and`(seq(`or`(F[i]>0, F[i]<0, F[i]=0), i=1..nops(F))))]; # The domain of res
solve(res) assuming `or`(seq(x in R[i],i=1..nops(R)));

The reason for the error (when calculating the derivative at a specific point) is a premature calculation. Here are 3 correct ways, starting with the best one:

restart;
y := t -> t^2: 
z := D(y);
z(1);

# or

y := t -> t^2:
z:=s->eval(diff(y(t),t), t=s);
z(1);

# or

y := t -> t^2:
z:=unapply(diff(y(t),t), t);
z(1);

restart;

S:=sum(1/x[i],i=1..5);
# or
S:=add(1/x[i],i=1..5);
# or
S:=`+`(seq(1/x[i],i=1..5));

That you ask in the additional question is called a polynomial decomposition (representation of a polynomial as a composition of two polynomials of lesser degree). This is (in general) impossible for arbitrary polynomials  g  and  h . If   f@g=h , then obviously the necessary condition  degree(h)=degree(f)*degree(g)  should be true that is not fulfilled in your example. Maple has the  compoly  command to decompose polynomials. See example below:

restart;
g:=unapply((x+1)^2, x);
f:=unapply(2*x^2+x+3, x);
h:=expand((f@g)(x));
compoly(h, x); 

expand(subs(%[2], %[1]));  # Check

In this example, we see that such a decomposition (if it is possible) is not unique.

Use the  tickmarks  option  for this. See help on  tickmarks  for details .
An example:

plot(x^2, x=-2..2, tickmarks=[0,0]);

Try spline fit. The approximation is so good that the 2 curves merge into one. Small differences can only be seen at narrower intervals:

restart;
Y := Vector([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 3, 8, 12, 12, 22, 30, 40, 44, 51, 65, 70, 97, 111, 131, 135, 174, 184, 210, 214, 232, 238, 254, 276, 285, 305, 318, 323, 343, 373, 407, 442, 493, 542, 627, 665, 782, 873, 981, 1095, 1182, 1273, 1337, 1532, 1728, 1932, 2170, 2388, 2558, 2802, 2950, 3145, 3526, 3912, 4151, 4399, 4641, 4787, 4971, 5162, 5445, 5621, 5959, 6175, 6401, 6677, 7016, 7261, 7526, 7839, 8068, 8344, 8733, 8915, 9302, 9855, 10162, 10819, 11166, 11516, 11844, 12233, 12486, 12801, 13464, 13873, 14554, 15181, 15682, 16085, 16658, 17148, 17735, 18480, 19147, 19808, 20244, 20919, 21371, 22420, 22614, 23298, 24077, 24567, 25133, 25694, 26484, 27110, 27564, 28162, 28711, 29268, 29789, 30249, 30784, 31323, 31987, 32558, 33153, 33616, 34259, 34854, 35454, 36107, 36663, 37225, 37801, 38344, 38948, 39539, 39977, 40532, 41180, 41804, 42208, 42689, 43151, 43537, 43841, 44129, 44433, 44890, 45244, 45687, 46140, 46577, 46867, 47290, 47743, 48116, 48445, 48770, 49068, 49485, 49895, 50488, 50964, 51304, 51905, 52227, 52584, 52800, 53021, 53317, 53477, 53727, 53865, 54008, 54247, 54463, 54588, 54743, 54905, 55005, 55160, 55456, 55632, 55829, 56017, 56177, 56256, 56388, 56478, 56604, 56735, 56956, 57145, 57243, 57437, 57613, 57724, 57849, 58062, 58198, 58324, 58460, 58647, 58848, 59001, 59127, 59287, 59345, 59465, 59583, 59738, 59841, 59992, 60103, 60266, 60430, 60655, 60834, 60982, 61194, 61307, 61440, 61558, 61630, 61667, 61805, 61882, 61930, 61992, 62111, 62224, 62371, 62521, 62691, 62853, 62964, 63036, 63173, 63328, 63508, 63731, 63790, 64090, 64184, 64336, 64516, 64728, 64884, 64996, 65148, 65305, 65693, 65693, 65839, 65982, 66228, 66230, 66439, 66607, 66805, 66974, 67220, 67330, 67412, 67557, 67838, 67960, 68303, 68627, 68937, 69225, 69645, 70195, 70609, 71344, 72140, 72757, 73175, 73394, 74132, 75062, 76207, 77013, 77933, 78434, 78790, 79789]):
N := numelems(Y): 
X := evalf(Vector([seq(1 .. N)])):
F:=CurveFitting:-Spline(X,Y, x):
P:=plot(F, x=0..N, color=red):
Q:=plot(X,Y, color=blue):
plots:-display(P,Q, size=[1000,350]);
plots:-display(P,Q, size=[1000,350], view=[150..160,40000..45000]);

Replace  D(y)(x)  by  diff(y(x), x) . Then the reason for the error becomes clearer:

restart;
eq:=exp(x)*sin(y(x))-3*x^2+(exp(x)*cos(y(x))+(1/3)/y(x)^(2/3))*(diff(y(x),x)) = 0;
solve(eq, y(x));
dsolve(eq, y(x), implicit);

restart;
v := a+2*H*Q/T; 
v := subs(H = T*x, v);

plot([subs(a = 0.1, Q = 0.5, v), subs(a = 0.2, Q = 0.4, v), subs(a = 0.3, Q = 0.35, v), subs(a = 0.4, Q = 0.3, v)], x = 1 .. 10, colour = [black, red, blue, green], axes = boxed);
 
plot([subs(a = 0.1, Q = 0.5, v), subs(a = 0.2, Q = 0.4, v)], x = 1 .. 10, style = pointline, colour = [black, red], symbol = [solidcircle, asterisk], symbolsize = 15, axes = boxed, numpoints = 20)

Try without the  parse  command:

restart;
P := plot(sin(x)-x+(x^3)/3!, x=-2*Pi..2*Pi) :
L := plots[textplot]( [ 3 , 17 , "sin(x)-x+(x^3)/3!"],  axes = none) :
plots[display](P ,L, axes = normal)

or

restart;
P := plot(sin(x)-x+(x^3)/3!, x=-2*Pi..2*Pi) :
L := plots[textplot]( [ 3 , 17 , sin(x)-x+(x^3)/`3`!],  axes = none) :
plots[display](P ,L, axes = normal)

This is the typical behavior of the  solve  command when solving transcendental equations. For symbolic (exact) representation of roots use the  identify  command:

restart;
Student:-Calculus1:-Roots(sin(x)-3*x/Pi, x=-Pi/4..Pi/4);
identify(%);

J4James There are many ways to smooth curves. Unfortunately your piecewise-curve (the red curve) has a discontinuity at the point x0=0.000588 (see the code below). Therefore, in any case, the approximation near point  x0  will not be very good. Below is my attempt of an approximation using the spline of degree 3 (the blue curve):

restart;

nuhf:=0.294*10^(-6):
Rehf:=Vhf*H/nuhf:
Vhf:=0.5:
rhohf:=965.31:
Lhs:=2.5*10^(-3):
DP1:=f1*Lhs*rhohf*Vhf^2/(2*H):
f1:=4*18.3/Rehf:
Whs:=1:
mhf:=rhohf*Vhf*Whs*H:
Rehf:=Vhf*H/nuhf:
PL1:=DP1*mhf/rhohf:
DP2:=f2*Lhs*rhohf*Vhf^2/(2*H):
f2:=0.3164/(Rehf^(0.25)):
PL2:=DP2*mhf/rhohf:
x0:=fsolve(Rehf=1000,H);
  
F:=unapply(piecewise(Rehf<1000, PL1, Rehf>=1000,PL2),H);
e:=5*10^(-5):
data:=[seq([H,F(H)],H=0.0001..0.0005,0.0001),[x0-e,F(x0-e)],[x0+e,F(x0+e)], seq([H,F(H)],H=0.0007..0.002,0.0001)];

P:=CurveFitting:-Spline(data, H);

A:=plot(F, 0.0001..0.002, color=red, discont):
B:=plot(P, H=0.0001..0.002, color=blue):

plots:-display(A,B, size=[500,500]);

Use  piecewise  instead of  if ... fi  :

restart;

nuhf:=0.294*10^(-6):
Rehf:=Vhf*H/nuhf:
Vhf:=0.5:
rhohf:=965.31:
Lhs:=2.5*10^(-3):
DP1:=f1*Lhs*rhohf*Vhf^2/(2*H):
f1:=4*18.3/Rehf:
Whs:=1:
mhf:=rhohf*Vhf*Whs*H:
Rehf:=Vhf*H/nuhf:
PL1:=DP1*mhf/rhohf:
DP2:=f2*Lhs*rhohf*Vhf^2/(2*H):
f2:=0.3164/(Rehf^(0.25)):
PL2:=DP2*mhf/rhohf:

procd:=proc(H) 
piecewise(Rehf<1000, PL1, Rehf>1000, PL2);
end:

plot(procd(H),H=0..0.002);