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

Since there are a large number of equivalent forms (as vv already mentioned), the simplest way is probably a completely manual solution or a Maple solution after manually splitting this expression into 3 parts:

B:=2*x^3 - 6*x*y*z + 2*y^3 + 2*z^3;
C:=factor(B);
op(1,C)*op(3,C);
``(x^2-2*x*y+y^2)+``(x^2-2*x*z+z^2)+``(z^2-2*z*y+y^2);
factor~(%);
expand(%);

Do not use the same name for the variable of integration and the limits of integration:

Download Appr_Int.mw

 Aragorn 15  I have corrected something in your worksheet. The procedure  orbit  now works successfully for  mode=1  and  mode=2  . 
 

Download Project_(1)_new1.mw

restart;
plots:-arrow([1,2], scaling=constrained, caption="first:1 \n second:2", captionfont=[times,16]);

What does  f  and  g  mean in your code? These are undefined. You probably meant that  f = Fresnelf  and  g= Fresnelg . 

Here's the corrected code:

restart;
a := 0;  
b := 1;
f:=Fresnelf:
g:=Fresnelg: 
p := plot({f(x), g(x)}, x = a .. b);
A := [a, g(a)];
B := [a, f(a)];                    
C := [b, g(b)];
D_ := [b, f(b)];                  
left_line := plot([A, B]);
right_line := plot([D_, C]);
plots[display](p, left_line, right_line, view = [0 .. b + 0.5, 0 .. 0.5]);

The left_line is just one point, because  A=B .

Alternatively, you can use the  alias  command as a shorthand:

alias(f=Fresnelf):
alias(g=Fresnelg): 

See help on  alias  command for details .

This will be a curve in space obtained from the intersection of the inclined plane  F  and the cylindrical surface  G:

restart;
F:=z=x+y;
G:= x^5+x^4*y+x^3*y+2*x^2*y^2+x*y^5+y^5-1=0;
A:=plots:-intersectplot(F, G, x=-5..5, y=-5..5, z=-5..5, thickness=3, numpoints=64000):
B:=plot3d(x+y, x=-5..5, y=-5..5, style=surface, color="LightBlue", view=-5..5):
plots:-display(A, B, axes=normal, orientation=[-40,60], lightmodel=light2);  

restart;
ls:=["a", "b", "c", "d", "e"]:
selectremove(t->ListTools:-Search(t,ls)::odd, ls);

The easiest way to solve the problem is to show your graphs in the  x  and  y ranges (using the  view  option), where the differences are clearly visible:

restart:
with(plots):
f1:=x->0.95*x;
f2:=x->0.99*x;
f3:=x->0.991*x;
p1:=plot(f1(x),x=0..10,color=red):
p2:=plot(f2(x),x=0..10,color=blue):
p3:=plot(f3(x),x=0..10,color=green):
P1:=display(p1,p2,p3, caption="Plot 1", captionfont=[times,16]); 
P2:=display(p1,p2,p3, view=[9.5..10,9.2..9.7],caption="Plot 2", captionfont=[times,16]);
display(<P1 | P2>);

Download plotting1.mw

Edit.

Maple finds the answer through elliptic functions. Unfortunately, the result is quite cumbersome:

restart;
y:=Int((k/2*x^2+x+c)/sqrt(lambda^2*(1+k*x)^4 - (k/2*x^2+x+c)^2), x);
value(y):
simplify(%);

Download Int.mw

Edit. In my opinion, the value of this overall result is rather dubious. With some parameters and variable  x  values, you can get  0 or negative numbers under the root. It is much more reliable to first set the values of the parameters, and only then calculate the integral.