This is a continuation of your last question. You should have given a link to that.
But here it is:
https://mapleprimes.com/questions/222957-Need-Help-Debbuging-My-Script-Howto-Get-My-Cdv

restart;
#mu := 0.18e-4; rho := 1.2; 
d := .2; Co := .4; # You may want to wait assigning to mu and rho
FDx := (3*Pi*mu*d+(1/8)*Pi*Co*rho*d^2*v)*vx;
FDy := (3*Pi*mu*d+(1/8)*Pi*Co*rho*d^2*v)*vy;
FD:=(3*Pi*mu*d+(1/8)*Pi*Co*rho*d^2*v)*v; #Using that v = sqrt(vx^2+vy^2)
FDD := (1/8)*Cdv*Pi*d^2*v^2; 
solve(FD=FDD,Cdv);
Cd:=unapply(expand(%),v);
## Testing
Cd(5);
mu := 0.18e-4; rho := 1.2;
plot(Cd(v),v=0..1,Cd=0..2);
############
Cd(v);
solve(R=v*d*rho/mu,{v});
CdR:=eval(Cd(v),%);
plot(CdR,R=1..10000);
plots:-loglogplot(CdR,R=1..10000);

The double logarithmic plot:

 

I would strongly recommend that in the future you delete all output before saving your document.
Use
Edit/Remove Output/From Worksheet.
If you have to hand in your homework as an mw-file with output, then I suggest that you save that version under a different name than your working copy.

Here is a solution found by the method of characteristics.
I'm following Fritz John, Partial Differential Equations, 1974, pp. 9-15.
I'm not initially finding a solution satisfying u(x,1/x) = 0, though.
At the end I give that problem a try.
 

restart;
PDE := (y-u(x, y))*diff(u(x, y), x)+(u(x, y)-x)*diff(u(x, y), y) = x-y;
#Following Fritz John, Partial Differential Equations, 1974: pp. 9-15.
a:=(x,y,z)->y-z;
b:=(x,y,z)->z-x;
c:=(x,y,z)->x-y;
sys:={diff(x(t),t)=a(x(t),y(t),z(t)),diff(y(t),t)=b(x(t),y(t),z(t)),diff(z(t),t)=c(x(t),y(t),z(t))};
dsolve(sys); #Same as the HINT=strip result
init:=[x0(s),y0(s),z0(s)];
#Condition (11) in Fritz John p. 11:
cond11:=diff(y0(s),s)*a(op(init))-diff(x0(s),s)*b(op(init))<>0;
#Choosing:
init1:={x0(s)=s,y0(s)=s+1,z0(s)=s};
eval(cond11,init1); #OK
# So
ics:={x(0)=s,y(0)=s+1,z(0)=s};
char:=dsolve(sys union ics);
XYZ:=subs(x(t)=x,y(t)=y,z(t)=z,char);
Z,XY:=selectremove(has,XYZ,z);
solve(XY,{s,t},explicit);
## The solution we find is then this:
U:=simplify(subs(%,rhs(op(Z))));
pdetest(u(x,y)=U,PDE); # 0
## The polynomial RD determines the domain of definition:
RD:=-9*x^2+18*x*y-9*y^2+12;
Student:-Precalculus:-CompleteSquare(RD,x);
plots:-implicitplot(RD,x=-3..3,y=-3..3);
plot3d(U,x=-2..2,y=x-2/sqrt(3)..x+2/sqrt(3));

The solution found above can be written  u(x,y) = ( x+y-sqrt(4-3*(x-y)^2))/2.

Now trying the problem with u(x,1/x) = 0. The result U2 is huge.
 

## Proceeding as in the simpler case until this point:
init2:={x0(s)=s,y0(s)=1/s,z0(s)=0}; #changed
eval(cond11,init2);
## New choice of ics:
ics2:={x(0)=s,y(0)=1/s,z(0)=0};
char2:=dsolve(sys union ics2);
XYZ:=subs(x(t)=x,y(t)=y,z(t)=z,char2);
Z,XY:=selectremove(has,XYZ,z);
solve(XY,{s,t},explicit): #Huge
U2:=simplify(subs(%,rhs(op(Z))),size): #Huge
#pdetest(u(x,y)=U2,PDE); # Does it finish?
eval(U2,y=1/x):
simplify(%) assuming x>0,x<1; # 0 so u(x,1/x) = 0.
plot3d(U2,x=0.001..1,y=1.001..3);
length(U2); # 40780

The expression u is positive, so abs(u) = u:
 

u:=-((-(1/20)*sqrt(5)-1/20)*cos((1/5)*Pi)-(1/20)*sqrt(2)*sqrt(5-sqrt(5))*sin((1/5)*Pi))*24^(1/5)*5^(3/5)-I*((1/20)*sqrt(2)*sqrt(5-sqrt(5))*cos((1/5)*Pi)+(-(1/20)*sqrt(5)-1/20)*sin((1/5)*Pi))*24^(1/5)*5^(3/5);
is(u>0);

A more detailed version:
 

restart;
u:=-((-(1/20)*sqrt(5)-1/20)*cos((1/5)*Pi)-(1/20)*sqrt(2)*sqrt(5-sqrt(5))*sin((1/5)*Pi))*24^(1/5)*5^(3/5)-I*((1/20)*sqrt(2)*sqrt(5-sqrt(5))*cos((1/5)*Pi)+(-(1/20)*sqrt(5)-1/20)*sin((1/5)*Pi))*24^(1/5)*5^(3/5);
is(u>0);
abs(u);
evalf(u);
u1:=abs(u-1);
u1-(1-u);

I have no problem reading the Norwegian. I think that Cd should be a function of the speed v, not the velocity (vx,vy):
 

restart;
#mu := 0.18e-4; rho := 1.2; 
d := .2; Co := .4; # You may want to wait assigning to mu and rho
FDx := (3*Pi*mu*d+(1/8)*Pi*Co*rho*d^2*v)*vx;
FDy := (3*Pi*mu*d+(1/8)*Pi*Co*rho*d^2*v)*vy;
FD:=(3*Pi*mu*d+(1/8)*Pi*Co*rho*d^2*v)*v; #Using that v = sqrt(vx^2+vy^2)
FDD := (1/8)*Cdv*Pi*d^2*v^2; 
solve(FD=FDD,Cdv);
Cd:=unapply(expand(%),v);
## Testing
Cd(5);
mu := 0.18e-4; rho := 1.2;
plot(Cd(v),v=0..1,Cd=0..2);

Just try pdsolve without condition and see what you get (or don't get).
Here I'm wrapping the call in timelimit because of my limited patience:
timelimit(60,pdsolve(PDE));

You probably already observed that u(x,y) = -x-y solves the PDE, but not u(x,1/x) = 0.

solve(sin(x),x,AllSolutions);

_Z1 (or _Z2, _Z3, ...) represents an arbitrary integer.

You could have mentioned that a related (but different) question was asked by you in
https://mapleprimes.com/questions/222830-How-To-Solve-Rational-Exponent-Differential#answer241834

In that link Rouben Rostamian gave an exact solution.
I gave a numerical solution, but to a revised equation where y(x)^0.337 was replaced by abs(y(x))^0.337.
That numerical solution was negative, whereas Rouben's exact solution was nonnegative.
In the present situation my numerical approach on the similarly revised equation works as before and delivers a negative solution.
Although probably not of interest I give the code here:
 

restart;
Digits:=15:
ode:=diff(y(x),x,x)-0.00003019*abs(y(x))^0.337=9.542*10^(-13)*x; 
bcs:=y(0)=0,D(y)(2945)=0;
res:=dsolve({ode,bcs},numeric,method=bvp[midrich],approxsoln=[y(x)=tanh(x)],initmesh=2048);
plots:-odeplot(res,[x,y(x)]);

This problem reminds me strongly of a problem that was answered beatifully by Rouben Rostamian.
https://mapleprimes.com/questions/222830-How-To-Solve-Rational-Exponent-Differential#answer241834
He used the fact that the initial value problem for the ode considered there had nonuniqueness for y(0) = y0 for some y0.
The same can be done here. Your ode has y(x) = A as a constant solution. It has nonuniqueness if y(0) = A is required.
This will be used to piece a solution together that has a constant part in the middle.
The solution is valid for 0 < A <= 2 (and 0 <= x <= 5).

restart;
ode:=y(x)*sqrt(1+diff(y(x),x)^2)-y(x)*diff(y(x),x)^2/sqrt(1+diff(y(x),x)^2)=A;
ode0:=normal(ode);
odetest(y(x)=A,ode0);
res:=solve(ode0,{diff(y(x),x)});
sol1:=dsolve(res[2] union {y(0)=2});
xA1:=solve(rhs(sol1)=A,x);
sol2:=dsolve(res[1] union {y(5)=5});
xA2:=solve(rhs(sol2)=A,x);
y(x)=piecewise(x<xA1,rhs(sol1),x<xA2,A,rhs(sol2)) assuming A>0,A<2;
solp:=combine(expand(%));
plots:-animate(plot,[rhs(solp),x=0..5],A=0..2);

Comment. By using ode0 we see that A must be positive since y(0) = 2 > 0.
Having noticed that it follows from ode0 that y(x) >= A for all x.
Thus in order for a piecewise solution as above to exist it must begin by going down (not up).

Start with interface( imaginaryunit=j) or some other name different from your variables since I is the imaginary unit by default in Maple.
A multiplication sign was missing after tau in lambda.
 

restart;
interface(imaginaryunit=j);
lambda := k*tau*(C*Upsilon+I)/N; #* was missing after tau
eqn1 := (1-p)*Pi+phi*V+delta*R-(mu+lambda+vartheta)*S;
eqn2 := p*Pi+vartheta*S-(epsilon*lambda+mu+phi)*V;
eqn3 := rho*lambda*S+rho*epsilon*lambda*V+(1-q)*eta*I-(mu+beta+chi)*C;
eqn4 := (1-rho)*lambda*S+(1-rho)*epsilon*lambda*V+chi*C-(mu+alpha+eta)*I;
eqn5 := beta*C+q*eta*I-(mu+delta)*R;
Equilibria := solve({eqn1 = 0, eqn2 = 0, eqn3 = 0, eqn4 = 0, eqn5 = 0}, {C, I, R, S, V});

The result looks perfectly fine to me. Just to replace your image with code in text form:
 

restart;
f:=1/s^2+(exp(s)-2)/s/(exp(s)-1);
g:=inttrans[invlaplace](f,s,t);
plot(g,t=0..5);
inttrans[laplace](g,t,s); # back to f
?frac  ## Help for frac

To understand how the answer can involve frac, look at this:
 

## f can be written as the infinite series:
fs:=1/s^2+1/s-Sum(exp(-n*s)/s,n=1..infinity);
## For illustration try this
inttrans[invlaplace](exp(-n*s)/s,s,t) assuming n>0; 
## Now transform fs:
gSum:=inttrans[invlaplace](fs,s,t);
plot(g-gSum,t=0..10,thickness=3); # zero

I'm assuming that by pi you meant the usual mathematical constant. That is spelled Pi in Maple.
pi in Maple is just another name. So I changed pi to Pi.
With this change there are no undetermined parameters in the system, so you can only give two boundary conditions.
I left out the one involving the derivative.
 

restart;
mu:=0.0025;
A:=1*10^(-20); Aprim:=4.46*10^(-20);
eta:=1.003*10^(-3);
V:=2*10^(-9) ;
ode:=mu*h(x)*diff(h(x),x,x)+mu*A/(6*Pi*h(x)^2)-eta*V; # pi ---> Pi
ode1:=mu*h(x)^3*diff(h(x),x,x)-eta*V*h(x)^2+mu*A/(6*Pi)=0; # pi ---> Pi
ic1:=h(0)=3/1000000000, h(-0.20e-4)=1e-9;
#D(h)(-0.20e-4)=0 left out
dsol1 := dsolve({ic1, ode}, numeric,maxmesh=256);
plots:-odeplot(dsol1,[x,h(x)]);

For a nonpolynomial equation or for systems of more than one equation fsolve only finds one solution at a time.
Start by plotting.
 

eqn1:=-x*y^2+4*x=5;
eqn2:=(1/3)*x^3+y^2=1;
plots:-implicitplot([eqn1,eqn2],x=-10..10,y=-10..10,gridrefine=1);
sol1:=fsolve({eqn1,eqn2},{x=-10..10,y=-10..10}); # As you did.
## Use the plot to find approximate solutions for a guess in the following:
sol2:=fsolve({eqn1,eqn2},{x=1,y=0.4});
sol3:=fsolve({eqn1,eqn2},{x=-3,y=3});
sol4:=fsolve({eqn1,eqn2},{x=-3,y=-3});

Since I is the imaginary unit in Maple you cannot assign to I(a).
But you can just use a name like A:
 

A:=int(exp(-x^2)*sin(a*x),x=0..infinity);
plot(A,a=-10..10);

If you don't like the appearance of the imaginary unit in the answer A, then you can convert it to 'dawson':
 

convert(A,dawson);
B:=combine(%);

The answer B is very simple when using dawson: B = dawson(a/2).
The definition of dawson is (use ?dawson):
dawson(x) = exp(-x^2)*int(exp(t^2), t = 0 .. x).

Because of the possibilty that imaginary results may show up because of your exponent and when y(x) <0 I changed y(x) to abs(y(x)).

Also I asked for a midpoint method (here midrich).
Thirdly I didn't have any initial success with the range 0..2945 so just to get something, I changed it to 0..450:

ode:=diff(y(x),x,x)-0.00003019*abs(y(x))^0.337=0; 
#bcs:=y(0)=0,D(y)(2945)=0.0116; 
bcs:=y(0)=0,D(y)(450)=0.0116; ## Try this first
res:=dsolve({ode,bcs},numeric,method=bvp[midrich]);
plots:-odeplot(res,[x,y(x)]);

Trying continuation in the exponent and using the full interval: Success!
 

restart;
ode:=diff(y(x),x,x)-0.00003019*abs(y(x))^0.337=0; 
bcs:=y(0)=0,D(y)(2945)=0.0116;
ode2:=evalindets(ode,`^`,x->op(1,x)^(c*op(2,x)));
res:=dsolve({ode2,bcs},numeric,method=bvp[midrich],continuation=c,maxmesh=8192);
plots:-odeplot(res,[x,y(x)]);
plots:-odeplot(res,[x,diff(y(x),x)]);