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

I don't know why  PDEtools:-Solve  doesn't return any result, but  explicit option helps:

PDEtools:-Solve(res, explicit);

Obviously, both solutions  sol  and  sol2  are the same for any number  x<>0 (that is, these are not exactly the same). For  x=0  the first form  sol  does not provide an obvious solution  y=-1. In this sense, the second form is better of course. But this is the design and so all the questions to the developers.

restart;
eq:=x*sin(x*y)+y+1=0;
sol:=solve(eq,y);
eval(%, x=0);
                  eq := x sin(x y) + y + 1 = 0
                            /      2            \
                      RootOf\_Z + x  sin(_Z) + x/
               sol := ---------------------------
                                   x             
Error, numeric exception: division by zero

sol2:=RootOf(x*sin(_Z*x)+_Z+1);
eval(%, x=0);
              sol2 := RootOf(x sin(_Z x) + _Z + 1)
                               -1
 

It seems that this function has discontinuities (or is it rounding errors)?

restart;

f := x->int((sin(t)^2)^exp(t), t = x .. infinity, digits=15, method = _d01amc, numeric);
plot(f, -5 .. 5, thickness=2, color=red);
plot(f, 3 .. 3.3, thickness=2, color=red);

Edit. Of course, there should be no discontinuities. The method below gives a smooth plot:

restart;

f := x->int((sin(t)^2)^exp(t), t = x..5,  numeric) + int((sin(t)^2)^exp(t), t = 5 .. infinity, method = _d01amc,  numeric);
plot(f, -5 .. 5, thickness=2, color=red);
plot(f, 3 .. 3.3, thickness=2, color=red);

To get a better look at this outlandish solid - here is its animation:

restart;
alpha[1] := sin(t):
alpha[2] := -1/2 - cos(t):
alpha[3] := 0:
beta[1] := 0:
beta[2] := 1/2 - cos(t)/(1+cos(t)):
beta[3] := sqrt(1 + 2*cos(t))/(1+cos(t)):
x := (1-m)*alpha[1] + m*beta[1];
y := (1-m)*alpha[2] + m*beta[2];
z := (1-m)*alpha[3] + m*beta[3];
P:=plot3d([[x,y,z], [x,y,-z], [z,-y,x], [-z,-y,x]], m=0..1, t=-Pi/2..Pi/2, style=surface, scaling=constrained, color=[red,green,blue,yellow]):

F:=t->plottools:-rotate(P,t,[[0,0,0],[0,0,1]]):
plots:-animate(F, [t], t=0..2*Pi, frames=90, axes=normal, view=[-2..2,-2..2,-1.5..1.5]);