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.

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