I think the bad news is that there will necessarily be a singularity in the interior when using conditions at x = Pi.
Thus since int(r(y),y=0..x) appears in the original formulation the conclusion is that the problem has no solution.
In the following I use the revised values you gave in your second worksheet.
First the code:
 

restart;
Digits:=15:
with(plots);
R0 := 10^(-5);
C := 1;
m := 1/100;
u1 := 1000;
sys := R(x) = 2*(-r(x)*diff(r(x), x)*cos(x)+r(x)^2*sin(x)+diff(r(x), x)^2*sin(x)-sin(x)*r(x)*diff(r(x),x,x))/(r(x)^4*sin(x)), 
       diff(R(x),x,x)+cot(x)*diff(R(x), x) = -(1/2)*r(x)^2*(R0^2-R(x)^2)-r(x)^2*C^2*m^2*exp(-2*m*A(x))/(2*u1), 
       diff(A(x), x) = r(x);
## I prefer this version:
SYS:=solve({sys},{diff(r(x),x,x),diff(R(x),x,x),diff(A(x),x)});
## I set the condition on A(x) at x = Pi, but of course we don't know its value.
## We let it be a parameter here called APi.
## 
condA := R(Pi) = 1/500^2, r(Pi) = 500, D(R)(Pi) = 0, D(r)(Pi) = 0, A(Pi) = APi;
## Now use the parameters option in dsolve/numeric:
res:=dsolve(SYS union {condA},numeric,parameters=[APi],output=listprocedure);
## Extracting the procedure for finding A(x) for given x:
Aproc:=subs(res,A(x));
## Setting as an example the parameter APi to 10^4:
res(parameters=[10^4]);
## The plots:
plots:-odeplot(res,[x,r(x)],0..Pi);
plots:-odeplot(res,[x,R(x)],0..Pi);
plots:-odeplot(res,[x,A(x)],0..Pi);
## To explore further we use this procedure, which takes APi and x values as input.
## If out=Aproc(xx) is defined for the given APi=aa then that  is the output, otherwise the output is -10^3.
## The most important is what is printed as messages:
##
Q:=proc(aa,xx) local out;
   if not type([aa,xx],list(realcons)) then return 'procname(_passed)' end if;
   res(parameters=[aa]); 
   try 
     out:=Aproc(xx);
     printf("Ok at APi = %f and x = %f\n",aa,xx);
   catch:
     printf(cat(StringTools:-FormatMessage( lastexception[ 2 .. -1 ]),"  APi=%f and x = %f\n"),aa,xx);
     out:=-10^3
   end try;
   out
end proc;
##
Q(2000,1); #Just a single call to Q with APi = 2000 and x = 1.
## Plotting Q(2000,x) (i.e. APi = 2000 for all the x values):
plot(Q(2000,x),x=0..Pi,adaptive=false,numpoints=50);
##
## Now plotting Q(a,1) as a function of a:
plot(Q(a,1),a=1..10^4,adaptive=false,numpoints=50);
plot(Q(a,0.9),a=1..10^4,adaptive=false,numpoints=50);
## It seems to be clear that for APi > 9000 we have reached the maximal interval of definition:
## That interval appears to be 0.99622186 .. Pi.
## It doesn't appear to help to raise APi above 9000.
## The reason is that the exp-term containing A(x) contributes almost nothing for larger values of APi.
##Compare to what happens if we remove that exp-term and therefore also the ode for A(x) as we don't need it:
#Removing the A-part totally from SYS:
SYS0:=remove(has,eval(SYS,exp=0),A(x));
res0:=dsolve(SYS0 union {R(Pi) = 1/500^2, r(Pi) = 500, D(R)(Pi) = 0, D(r)(Pi) = 0},numeric);
odeplot(res0,[x,r(x)],0..Pi); #Notice the value 0.99622186
odeplot(res0,[x,R(x)],0..Pi);

Since you have the integral Ar(x) = int(r(x1), x1 = 0 .. x) you can introduce an ode for Ar:
diff(Ar(x),x) = r(x) with Ar(0) = 0.
You may have a serious singularity at least at r = 0 as will be exhibited in the code below:

restart;
R0 := 10^(-5): C := 1: m := 1: u1 := 1:
sys := R(x) = 2*(-r(x)*diff(r(x), x)*cos(x)+r(x)^2*sin(x)+(diff(r(x), x))^2*sin(x)-sin(x)*r(x)*diff(r(x), x$2))/(r(x)^4*sin(x)), 
      diff(R(x), x$2)+cot(x)*diff(R(x), x) = -(1/2)*r(x)^2*(R0^2-R(x)^2)-r(x)^2*T(x)/(2*u1);
T := x->C^2*m^2*exp(-2*m*Ar(x));
odeAr:=diff(Ar(x),x)=r(x); # Adding this ode
condAr:=Ar(0)=0; # Boundary value for Ar.
cond := R(Pi) = 1/500^2, r(Pi) = 500, D(R)(Pi) = 0, D(r)(Pi) = 0; #Your given boundary values.
##
#Solving the system sys, odeAr for the highest derivatives to see the singularity problem:
SYS:=solve({sys,odeAr},{diff(r(x),x,x),diff(R(x),x,x),diff(Ar(x),x)});
expand(SYS); #Possible problem at x=0 (and perhaps x=Pi)
F:=dsolve(SYS union {cond,condAr},numeric,method=bvp[midrich],initmesh=8192);

I get the error message:

Error, (in dsolve/numeric/bvp) powering may produce overflow

The ode for r(x) in the isolated version is
 

diff(r(x), x, x) = -(1/2)*r(x)^3*R(x)+r(x)-diff(r(x), x)*cos(x)/sin(x)+diff(r(x), x)^2/r(x);

You notice that the last two terms have denominators sin(x) and r(x), respectively. But you need to look at the numerators too, of course.
You could try to avoid initial problems by using an approximate solution where the initial guess for r(x) has r(0)<>0, D(r)(0) = 0, and satisfies the boundary condition D(r)(Pi) = 0 as required, but I haven't had any success with my first attempts in that direction.

 

Since the plot actually shows that the function a[1] is constant and equal to 0.1 you could plot the difference a[1](t)-0.1 as in this example where I use the logistic equation:
 

restart;
ode:=diff(x(t),t)=x(t)*(0.1-x(t));
res:=dsolve({ode,x(0)=0.1},numeric);
plots:-odeplot(res,[t,x(t)],0..10);
plots:-odeplot(res,[t,x(t)-0.1],0..10,caption="The difference between x(t) and 0.1");

The last plot:

I don't know exactly what was done in 3(d), but I would do the following, where the important part is to draw the orbit after it has had time to get close to the apparent attracting limit cycle. Thus I plot from t=200 to t=230:
 

restart;
with(plots):
EquationSet:= {diff(x(t),t) = -y(t) - z(t), diff(y(t),t) = x(t) + y(t)/5, diff(z(t),t) = 1/5 + (x(t) - 5/2)*z(t)};
Numsol1:=dsolve({EquationSet[],x(0)=0,y(0)=0,z(0)=0},[x(t),y(t),z(t)],numeric):
odeplot(Numsol1,[x(t),y(t),z(t)],200..230,numpoints=10000);
odeplot(Numsol1,[x(t),y(t)],200..230,numpoints=10000);
EquationSet2:= {diff(x(t),t) = -y(t) - z(t), diff(y(t),t) = x(t) + y(t)/5, diff(z(t),t) = 1/5 + (x(t) - 3)*z(t)};
Numsol2:=dsolve({EquationSet2[],x(0)=0,y(0)=0,z(0)=0},[x(t),y(t),z(t)],numeric):
odeplot(Numsol2,[x(t),y(t),z(t)],200..230,numpoints=10000);
odeplot(Numsol2,[x(t),y(t)],200..230,numpoints=10000);
## To produce an animation of the change from 2 to 3 via the value 5/2 I use the parameters option in dsolve/numeric.
##
EquationSetA:= {diff(x(t),t) = -y(t) - z(t), diff(y(t),t) = x(t) + y(t)/5, diff(z(t),t) = 1/5 + (x(t) - A)*z(t)};
res:=dsolve({EquationSetA[],x(0)=0,y(0)=0,z(0)=0},[x(t),y(t),z(t)],numeric,parameters=[A]);
Q:=proc(A) if not A::realcons then return 'procname(_passed)' end if;
      res(parameters=[A]);
      res
end proc;
animate(odeplot,[Q(A),[x(t),y(t),z(t)],200..230,numpoints=10000],A=2..3,frames=50,orientation=[0,20,30]);

You have Duffing's equation written as a first order autonomous system.
Instead of presenting the code (or in addition to that) you should explain what you are trying to do.
In the meantime you can have a look at this code, which uses the original second order nonautonomous ode:

restart;
ode := diff(x(t), t,t) = -2*gamma1*diff(x(t),t)-alpha*x(t)-beta*x(t)^3+F*cos(t);
ics:={x(0)=x0,D(x)(0)=x1};
params:=[alpha,beta,gamma1,F,x0,x1];
res:=dsolve({ode} union ics,numeric,parameters=params);
res(parameters= (params=~[-1,1,0.125,0.3,0,0])); #Setting parameters
plots:-odeplot(res,[t,x(t)],0..200,size=[1800,default],numpoints=1000);
plots:-odeplot(res,[x(t),diff(x(t),t)],0..200,numpoints=10000,frames=100);
plots:-odeplot(res,[x(t),diff(x(t),t),t],0..200,numpoints=10000,frames=100);

Above you see the first animation.

Addendum.
I guess you are trying to do something like the following, where I give the full code by making some changes to the code above:
 

restart;
Digits:=15:
ode := diff(x(t), t,t) = -2*gamma1*diff(x(t),t)-alpha*x(t)-beta*x(t)^3+F*cos(t);
ics:={x(0)=x0,D(x)(0)=x1};
params:=[alpha,beta,gamma1,F,x0,x1];
res:=dsolve({ode} union ics,numeric,parameters=params,output=listprocedure,maxfun=10^7,abserr=1e-14,relerr=1e-12);
X,X1:=op(subs(res,[x(t),diff(x(t),t)]));
res(parameters= (params=~[-1,1,0.125,0.3,0,0]));
numpts:=1000: #You had 50.
pts := seq([X(2*Pi*i+200), X1(2*Pi*i+200)], i = 1 .. numpts):
plot([pts],style=point,color=blue,scaling=constrained);

The plot resembles the one in Guckenheimer and Holmes, p. 87, fig. 2.2.4 (b).

You may experiment with this:
 

dsol5 := dsolve(  dsys3,numeric,stiff=true,maxfun=10^7); 
plots:-odeplot(dsol5, [t, r[1](t)], 0 .. 0.9e-3);

In your Maple code you are using programmer indexing as opposed to math indexing.

If you change to the latter you get:
 

restart:
with(CodeGeneration):

QLoc:=proc()
local Q:
Q:= Matrix(2,2):
Q[1,1] := 1E5:
Q[2,2] := 1E4:
Q[1,2] := 1E3:
Q[2,1] := 1E3:
return Q:
end proc:
 
Python(QLoc);

import numpy

def QLoc ():
    Q = numpy.mat([[0,0],[0,0]])
    Q[0][0] = 0.1e6
    Q[1][1] = 0.1e5
    Q[0][1] = 0.1e4
    Q[1][0] = 0.1e4
    return(Q)

Note: I don't know Python, so I can't tell you if this is any better.
You can read about the two kinds of indexing in Maple in the help page:
?Indexing Arrays, Matrices, and Vectors

Maybe you have the wrong value for x(0) = x0. Also you need a value for the derivative at zero.
Your ode can be written as a system in the usual way (see below). It is autonomous and has 3 equilibrium points:
(0,0), (sqrt(2),0),  and (-sqrt(2),0) ( with your values of a and b: are they correct?).
(0,0) is asymptotically stable, the other two are saddle points, thus unstable.
A solution starting with x(0) = 2 and x'(0) = 0 will not be defined for all t = 0..100: It "reaches" infinity in a finite time.

restart;
ode:=diff(x(t),t,t)+k*diff(x(t),t)=D(V)(x(t));
V:=x->a*x^2/2+b*x^4/4;
params:={a=-2,b=1,k=0.1};
ode;
sys:={diff(x(t),t)=y(t), diff(y(t),t)=-k*y(t)+b*x(t)^3+a*x(t)};
solve({y=0,-k*y+b*x^3+a*x},{x,y},explicit) assuming a<0,b>0;
E:=map(subs,[%],[x,y]);
eval(E,params);
J:=unapply(VectorCalculus:-Jacobian([y,-k*y+b*x^3+a*x],[x,y]),x,y);
J(op(E[1]));
LinearAlgebra:-Eigenvalues(%);
eval(%,params);
J(op(E[2]));
LinearAlgebra:-Eigenvalues(%);
eval(%,params);
eval(LinearAlgebra:-Eigenvalues(J(op(E[3]))),params);
DEtools:-DEplot(eval(sys,params),[x(t),y(t)],t=-10..100,[[x(0)=1,y(0)=0],[x(0)=1.5,y(0)=0],[x(0)=-1.5,y(0)=0]],
   x=-3..3,y=-3..3,linecolor=blue,stepsize=0.1);
res:=dsolve({eval(ode,params),x(0)=2,D(x)(0)=0},numeric);
plots:-odeplot(res,[x(t),diff(x(t),t)],-1.2..1.2,view=[0..10,-10..10],numpoints=10000);

The DEplot orbits:


The dsolve/numeric orbit:

First of all you have phi appearing as just the name phi, but also as the function phi(eta).
Since an actual number (like 3.57) can act as a constant function, i.e.  3.57(eta) is just 3.57 for all eta, this confusion won't matter in your case. But replace in Eq2 phi(eta) by phi anyway.
Secondly, there is no need for using continuation, but use maxmesh = 1024.
Thirdly, you are asking for dsol[h][i](0), but eta=0 is outside of the interval. You need to replace that with an eta belonging to the interval. Maybe the left endpoint eta = a = 0.1.
So the crucial line will be

dsol[h][i] := dsolve(eval(dsys, d = 1), numeric, maxmesh = 1024);

Finally, since all three values of a are the same and all three values of phi are the same there is obviously no need for the double loop.
This would do:
 

res:=dsolve(eval(dsys,{phi=0.05,a=0.1,d=1}), numeric,maxmesh=1024);
## Plotting
odeplot(res,[eta,f(eta)]);
odeplot(res,[eta,theta(eta)]);

The plot of theta:

By x^(1/3) is in Maple meant the principal cube root. There are 3 cube roots.
The principal root is the one which has principal argument of least absolute value. In case of ties the one with positive argument is chosen.
With this understanding your equation doesn't have a root:
 

restart;
u:=4+x^(1/3);
solve(u=0);
eval(u,x=r*exp(I*t)); #Polar decomposition: r>=0, t>-Pi, t<=Pi.
w:=simplify(evalc(%)) assuming r>0,t>-Pi,t<=Pi;
## We need the imaginary and the real parts to be zero
solve(sin(t/3)=0,t,allsolutions);
## Must choose Z=0 for t to belong to (-Pi, Pi]
eval(w,t=0); # Never zero for r >= 0

restart;
eq1 := diff(f(eta),eta,eta,eta,eta)+(2/(eta+k))*diff(f(eta),eta,eta,eta)-(1/(eta+k)^2)*diff(f(eta),eta,eta)+(1/(eta+k)^3)*diff(f(eta),eta)+B*(k/(eta+k)*(F(eta)*diff(f(eta),eta,eta,eta)-diff(F(eta),eta)*diff(f(eta),eta,eta)-diff(F(eta),eta,eta)*diff(f(eta),eta)-diff(F(eta),eta,eta,eta)*f(eta))+(k/(eta+k)^2)*(F(eta)*diff(f(eta),eta,eta)-2*diff(F(eta),eta)*diff(f(eta),eta)+diff(F(eta),eta,eta)*f(eta))-(k/(eta+k)^3)*(F(eta)*diff(f(eta),eta)+diff(F(eta),eta)*f(eta))-(beta/(eta+k))*(diff(f(eta),eta)+(eta/2)*diff(f(eta),eta,eta))-(beta/2)*(3*diff(f(eta),eta,eta)+(eta)*diff(f(eta),eta,eta,eta))+(A/(eta+k))*diff(f(eta),eta)+(A)*diff(f(eta),eta,eta));
eq2 := diff(F(eta),eta,eta,eta,eta)+(2/(eta+k))*diff(F(eta),eta,eta,eta)-(1/(eta+k)^2)*diff(F(eta),eta,eta)+(1/(eta+k)^3)*diff(F(eta),eta)+B*((k/(eta+k))*(diff(F(eta),eta)*diff(F(eta),eta,eta)-F(eta)*diff(F(eta),eta,eta,eta))+(k/(eta+k)^2)*(diff(F(eta),eta)*diff(F(eta),eta)-diff(F(eta),eta,eta)*F(eta))-(k/(eta+k)^3)*(F(eta)*diff(F(eta),eta))-(beta/(eta+k))*(diff(F(eta),eta)+(eta/2)*diff(F(eta),eta,eta))-(beta/2)*(3*diff(F(eta),eta,eta)+(eta)*diff(F(eta),eta,eta,eta)));
bc:=F(0)=S,D(F)(0)=lambda,D(F)(N)=0,D(D(F))(N)=0,f(0)=0,D(f)(0)=0,D(D(f))(0)=1,D(f)(N)=0,D(D(f))(N)=0;
N:=6;lambda:=-1;beta:=-2; k:=200;S:=2;cphi:=0.2; rhos:=3970; rhof:=1115;
B:=(1-cphi)^(2.5)*(1-cphi+cphi*(rhos/rhof));
A1:=dsolve({eq1,eq2,bc},numeric,output=operator,initmesh=128,approxsoln=[A=3.46,F(eta)=2*exp(-eta)-0.1,f(eta)=0.2*tanh(eta)]);
A1(0); #Here you see the value of A: 3.45824367723784
plots:-odeplot(A1,[eta,F(eta)]); 
plots:-odeplot(A1,[eta,f(eta)]);

To get you started here is the code for solving the system numerically. I'm using the constants from your previous question, which seems to deal with the same thing:
https://mapleprimes.com/questions/223001-Error-In-DEplot3D

restart;
ode1:=diff(x(t),t)=0.1*x(t)-0.00008*x(t)*y(t);
ode2:=diff(y(t),t)=0.00004*x(t)*y(t)-0.04*y(t);
ics:={x(0) = 1, y(0) = 0.1};
res:=dsolve({ode1,ode2} union ics,numeric,range=0..250);
plots:-odeplot(res,[x(t),y(t)],0..250,thickness=3,refine=1);
plots:-odeplot(res,[x(t),y(t),t],0..250,thickness=3,refine=1);

As you will surely know, the orbits of all solutions starting in the open first quadrant are closed.

There is a serious syntax error in your system of odes:
You have cos(varphi))(t) where you ought to have cos(varphi(t) ).
This error is repeated throughout the expression and also with sin.
There may be more problems. Since, however, you are asking a related question here:
https://mapleprimes.com/questions/223002-How-Plot-Phase-Portrait-On-A-Cylinder
I shall not look into this further.

The singularity is exhibited very clearly graphically in the plots:

odeplot(F,[theta,diff(mu(theta),theta)],0..4);
## and a close up:
odeplot(F,[theta,diff(mu(theta),theta)],2.9..3.2,view=-1..0);

Isolating the second derivative from the second ode:

combine(solve(sys[2],{diff(mu(theta),theta,theta)}));

you get the ode (in a short notation):
mu'' = -cot(theta)*mu' +1-exp(2*mu)*R/2
Now assume that mu, mu', and R are continuous at Pi (or at least have finite limits from the left at Pi).
Assume further that mu'(theta) -> mu1 for theta->Pi from the left and that mu1<>0.
Then it follows from the ode that mu''(theta) -> signum(m1)*infinity as theta -> Pi from the left.
From the second graph of mu' it appears that mu'(theta) clearly stays negative for theta < Pi.
 

The help page says:

The select function selects the operands of the expression e which satisfy the Boolean-valued procedure f, and creates a new object of the same type as e.

The type in your case is function, in Maple meaning an unevaluated function call as e.g. sin(2) or f(x).
The only operand of f(x) is x: Try op(f(x));
That you can get f by doing op(0,f(x)) is convenient, but doesn't make f an operand of f(x) in the usual sense.
As Kitonum has already answered: Wrap D[1](xx[0]) in a list or a set.

Details:

restart; 
op(f(x)); 
op(0,f(x)); 
remove(has, f(x), x); #The result must be a function since f(x) is. 
type(f(),function); # true
select(has, f(x),x); #The result must be a function since f(x) is. 
selectremove(has,f(x), x); #The results must be functions since f(x) is.
## Now wrap in { } :
remove(has,{f(x)},x); # The result must be a set since {f(x)} is a set. 
select(has,{f(x)},x); # The result must be a set since {f(x)} is a set. 
selectremove(has,{f(x)},x); # The results must be sets since {f(x)} is a set.