I toke a quick glance at your worksheet.
First you define D13 as one system, then you define it as another system. That second system surely is different from the first. So why expect the same?
I don't see any limit cycle in the first case (and in fact not even in the second case).
Your first system:

Your second system:

@Carl Love I agree about the fourier sine series except for a factor. I think it should be 
4*Sum(....)= ln((1+sin(2*Pi*x))/(1-sin(2*Pi*x))).

How is Trajectoire defined (if at all)?

The answer for X>0 is clearly wrong (in Maple 17), even if you start by using a concrete a.
For positive a and in a range on which a*x^3+1 is positive INT is an increasing function of X whether or not X>0 or X<0.

INT:=Int(1/sqrt(a*x^3+1),x=0..X);

value(eval(INT,a=0.1)) assuming X>0;
plot(%,X=0..1); #Decreasing
#For X<0 OK:
value(eval(INT,a=0.1)) assuming X<0;
plot(%,X=-1..0); #Increasing

@bwarga Sorry. I might as well have said that inits was

inits:=[[x(0) = 0.1e-1, y(0) = .99], [x(0) = -.1, y(0) = -.9], [x(0) = 1.1, y(0) = 0], [x(0) = 0, y(0) = .2], [x(0) = 0, y(0) = .6], [x(0) = .6, y(0) = 0], [x(0) = .75, y(0) = 1], [x(0) = .1, y(0) = .1], [x(0) = .5, y(0) = 1.0], [x(0) = -.5, y(0) = 1], [x(0) = .5, y(0) = -1], [x(0) = -.5, y(0) = -1], [x(0) = -0.1e-1, y(0) = .99], [x(0) = 0.1e-1, y(0) = -.99], [x(0) = -0.1e-1, y(0) = -.99], [x(0) = .5, y(0) = -1], [x(0) = -.5, y(0) = -1], [x(0) = 0.1e-1, y(0) = .9]];

#Thus with that and

with(DEtools);
L := -1.576674; MU := 0; DE13 := {(D(x))(t) = x(t)*(1+4*x(t)*x(t)-y(t)*y(t))+MU*y(t)*(x(t)*x(t)-.43*y(t)*y(t)-L), (D(y))(t) = y(t)*(1+x(t)*x(t)-.5*y(t)*y(t))+MU*x(t)*(x(t)*x(t)-.43*y(t)*y(t)-L)};

#you get

DEplot(DE13, [x(t), y(t)], t = -5 .. 2, inits,x=-5..5,y=-2..2,
stepsize = 0.1e-1, title = "phaseplane 3 prime plot", linecolor = black, thickness = 1);

@Carl Love Having taken a second look, I found unevaluated integrals in Maple 15 and 16.The reason I didn't see them the first time was that I just executed the whole worksheet with no error error messages appearing and the two plots produced in the 3 Maple versions (15,16,17) I have available here.
The unevaluated integral used for f(tau) in Maple 16 (and 15) doesn't get in the way of finding
diff(f(tau),tau) explicitly
nor would it be a problem to plot f(tau) (even if still unevaluated):
plot(f(tau),tau=0..5);

I don't know what the questioner meant by "the code does not converge".

I don't get any such Z1 in Maple 17,16,15.

I haven't updated the Physics package. I get no such error in Maple 17.02.

@Stephan This seems to me to be a bug.
You can use isolate instead:

isolate(eq[1],sigma[F](xi));

To see that it is a bug you can try

restart;
p:=piecewise(x<1,67,0);
solve(p=a/b(0),a);
solve(p=a/b,a);
solve(p=a/exp(1),a);


I'm using Maple 17 as you are. I don't know if this is fixed in Maple 18.
The bug is NOT present in Maple 15 nor in Maple 16.

Just for the record I shall submit an SCR (bug report).

@J4James I think I finally got what you want:

restart;
ODE:=diff(u(y),y)=y*m*exp(Lambda*(y*m)^2);
bcs:=D(u)(0)=0,u(h)=-1;
dsolve({ODE,bcs});
uu:=rhs(%);
us:=simplify(F-int(uu,y=0..h)) assuming Lambda>0;
h:=1+phi*cos(2*Pi*x):
F:=A-1:
Lambda:=2*B*W^2:
A:=2:x:=1:B:=3:W:=2:
uu;
U1:=evalf(convert(us,dawson));
p:=proc(phi1) if not type(phi1,numeric) then return 'procname(_passed)' end if;
   fsolve(eval(U1,phi=phi1)=0,m)
end proc;
p(0);
p(0.2);
U:=phi1->evalf(eval(uu,{phi=phi1,m=p(phi1)}));
U(0);
plot([U(0),U(0.2)],y=-1..1,axes=box,color=[red,green]);
plots:-animate(plot,[U(phi),y=-1..1,axes=box],phi=0..0.5);

@J4James
restart;
ODE:=diff(u(y),y)=y*m*exp(Lambda*(y*m)^2);
bcs:=D(u)(0)=0,u(h)=-1;
dsolve({ODE,bcs});
uu:=rhs(%);
us:=simplify(F-int(uu,y=0..h)) assuming Lambda>0;
h:=1+phi*cos(2*Pi*x):
F:=A-1:
Lambda:=2*B*W^2:
A:=2:x:=1:B:=3:W:=2:
U:=evalf(convert(us,dawson));
p:=proc(phi1) fsolve(eval(U,phi=phi1)=0,m) end proc;
p(0);
p(0.2);
plot(p,0..0.2,adaptive=false,numpoints=5);

But I don't understand where your variable y is coming from in your "goal".

@mehdi jafari If you change the output in your procedure to

Array(1..Fx,1..Fy,1..Ft,(i,j,k)->u[i,j,k],datatype=float)

then you can use matrixplot as in

plots:-matrixplot(Matrix(M[..,..,5]),labels=[x,y,u],caption=typeset("Time fixed at t =%1",5));

And an animation in time:

S:=seq(plots:-matrixplot(Matrix(M[..,..,i]),labels=[x,y,u],caption=typeset("Time t =%1",i)),i=1..21):
plots:-display(S,insequence);

Notice that time is just represented by the index.

@J4James I don't think that there is any way of getting a solution for m in closed form.
You need values for all parameters to be able to find the (unique) solution there is for m.

@Stephan Yes you can easily remove them using evalindets. But how did they get there in the first place?
It might be better to correct the problem at that stage.
However, here is the removal by evalindets:

evalindets(FuncA,list,op);

@Stephan 

Funct:=piecewise(xi*L<=L/2,x^2,L/2>xi*L,x);
Funct2:=eval(Funct,L=1); #The extremely simple solution for your specific example
evalindets(Funct,relation,s->map(x->x/L,s)); #the more complicated
evalindets(Funct,`<=`,s->map(x->x/L,s)); #restricted to type `<=`

A slightly more complicated example shows the broader scope of using evalindets:

G:=piecewise(xi*L<=L/2,x^2+L*x,L/2>xi*L,x+sin(x*L));
evalindets(G,relation,s->map(x->x/L,s));