Your worksheet only contains this input:
 

w := 1/2*(beta*eta^2*y + 2*eta^2*y + 2*A*beta + 2*A)/kx + C1*exp(k*y) + C2*exp(-k*y) + C3*exp(-I*k*y) + C4*exp(k*y*I);

I don't see any ode or BCs.

@mmcdara In order to avoid the existence of several local variants of a variable place all assumptions at the beginning, at least before making assignments involving the variables in the assumptions.
As in

assume(r > 0,R>0,R>=r,x__Q>0,d>-R-2*r,d<R);

But I much prefer using assuming rather than assume. Then you don't have that kind of problem.
The assume facility is not perfect either, of course.
Using assuming in your second attempt:
 

restart:
interface(rtablesize=10):
with(geometry):
_EnvHorizontalName := x: _EnvVerticalName := y:
point(P,0,0);
circle(A,[P, R]) assuming R>0;

@Zeineb If the result of the solve command is called res (say), then res will be a sequence having (as it turns out) 4 elements.
To see that just do:  nops( [res] );
res[1] and res[2] are the trivial solutions:
The next two both involve RootOf. (see ?RootOf).
You could say they are just reformulations of the two equations in your input.
You probably won't be able to simplify those in any significant way.
I'm afraid that your equations can only be solved numerically or in special cases where some of the parameters are known.
For a numerical solution you need to know all the parameters.
For the numerical solution you can use fsolve, but can also use the RootOf results as illustrated here:

params:=indets({firstequation,secondequation},name) minus {rho1,rho2,N}; # N treated separately
## Example: All parameters except N are 1 and N=2:
eval(res[3],(params=~1) union {N=2});
evalf(%); # result: {rho1 = .6288883452, rho2 = .6288883452}
eval(res[4],(params=~1) union {N=2});
evalf[20](%); # {rho1 = -.42264973081037423549+1.5788764000904431175*10^(-20)*I, rho2 = 1.}

@acer I couldn't find the option 'real' documented. Is it?

@ Sorry, I overlooked that.
Now I see that you have a boundary value problem on the interval 0 .. H-He.
The parameters option in dsolve is only available for initial value problems so the way you are handling the problem seems fine to me. The whole thing runs if you make the change to SYS as discussed above:
 

## All your code above the next line assumed here !
unassign(sv,sh,L,Jir,u):

sys_ode:=
 diff(sv(x),x)=rho0*g/L(x),
 diff(sh(x),x)=Fhp1*diff(L(x),x)/L(x)+Fhp2*diff(Jir(x),x),
 Fvp1*diff(L(x),x)+Fvp2*L(x)*diff(Jir(x),x)=rho0*g,
 diff(Jir(x),x)/Jir(x)=-Gp1*diff(L(x),x)/L(x)-Gp2*diff(Jir(x),x),
 diff(u(x),x)=-Mp*g/(Fvp1+Fvp2*L(x)*diff(Jir(x),x)/diff(L(x),x)):
##
dvars:=indets({sys_ode},specfunc(anything,diff));
SYS:=solve({sys_ode},dvars):
##
dt:=(Ts-Te)/nnp:
t:=Te:
H:=He:
lx0a:=Le:
Jirx0a:=Jire:
##
for i from 1 to nnp+1 by 1 do
   lprint(i,H-He); # Just to see H-He
   ics:=sv(H-He)=sve,sh(H-He)=she,L(H-He)=Le,Jir(H-He)=Jire,u(0)=0:
   sol:=dsolve(SYS union {ics},numeric,output=listprocedure):
   SV,SH,LA,JIR,U:=op(subs(sol,[sv(x),sh(x),L(x),Jir(x),u(x)])):
   tmy:=t/(10^6*365*24*60*60):
   lx0:=LA(0):
   Jirx0:=JIR(0):
   phix0:=1-(1-phi0)/Jirx0:
   shx0:=SH(0)/10^6:
   svx0:=SV(0)/10^6:
   pcx0:=(pc0*(ln(phix0)/ln(phi0))^m)/10^6:
   pvx0:=(pv0*(ln(phix0)/ln(phi0))^n)/10^6:

   up:=U(H-He):
   uh:=ue+up:
   fvpx0:=evalf(A*shx0+B*svx0-pvx0):

   if fvpx0>=0 then
     break
   else
     dH:=dt*(Mp/rho0+uh):
     H:=H+dH:
     t:=t+dt:
   end if:
end do:

This is obviously a bug. Not many, if any,  would have any use for that print.

I downloaded your worksheet. Ran it, and added an extra line.
Will try to upload with display.

Maple Worksheet - Error

Failed to load the worksheet /maplenet/convert/fibon2019_pa.mw .
 

Download fibon2019_pa.mw

@mstouby It will be, I'm sure.

Inserting the matrix on a 1D-input line works. 1D input is also called Maple Notation.
In my view 1D input is way superior to 2D.
This was another example.

@Nia Dutta How did you produce this "approximation"? It is not good anywhere except at one point, basically.

@rlopez Yes, I agree. But if sin(x)/x is the exact solution then that approximation is not only bad, but terrible:
 

numapprox:-infnorm(x^2+423*x^3/3+6789*x^4/4-sin(x)/x,x=0..2);

Value: 28287.54535

@Haqphy12015 Since your a2F(w) values are negative in the beginning the integral will necessarily be negative in the beginning. What to do about that must depend on the physical meaning of a2F(w) (assuming there is one).
Obviously using abs(a2F(w)) instead of a2F(w) will make the integral nonnegative and increasing, but so would e.g. a2F(w)^2.
Which physical meaning do w and a2F(w) have if any?

@tomleslie You are right. A person asking the same question twice (although somewhat differently phrased) at least ought to give us a link to the earlier question and tell us why he wasn't satisfied with the first answer.

@Carl Love Thanks. ApproximateInt fortunately doesn't shy away from symbolic input so it can be used to show your point:
 

restart;
S1:=Student:-Calculus1:-ApproximateInt(f(x),x=a..b,method=right,partition=1);
S2:=Student:-Calculus1:-ApproximateInt(f(x),x=a..b,method=left,partition=1);
factor((S1+S2)/2);
Student:-Calculus1:-ApproximateInt(f(x),x=a..b,method=trapezoid,partition=1);

@acer Thanks for reminding me of the package Interpolation, which is new to Maple 2018.
Also for pointing out that the procedural version of numerical integration is much faster than the nonprocedural that I was using.