@Bendesarts If the limit cycle is asymptotically stable, which it appears to be, then in both of your examples you should just plot the solution on an interval of short length T:  tmax-T..tmax.

(Your tmax values are way over done. You can use considerably smaller values for tmax. Try tmax/20 in both cases, e.g.).

In your first example:

res:=dsolve({EqSys[],ic[]},numeric,range=0..tmax);
plots:-odeplot(res,[x(t),z(t)],0..tmax,refine=1);
plots:-odeplot(res,[x(t),z(t)],tmax-6.3..tmax,refine=1,scaling=constrained);



In your second example:

plots:-odeplot(res,[u[1](t),v[1](t)],tmax-0.6..tmax, scaling = constrained);

@Cata For clarity let us not use the variable pi since it prints the same as Pi.
So consider the equation

eq:=arctan(x)+arctan(x/2)=q;
res:=solve(eq,x);
# You get a sequence of 2 solutions. This doesn't necessarily mean that there are two solutions.
Here the interpretation must be that one expression is used for certain q-values, the other for the rest of them.

Try plotting the left hand side of the equation (i.e. arctan(x)+arctan(x/2) ):
plot(lhs(eq),x=-25..25);
#Observe also that the range of arctan(x)+arctan(x/2) is the open interval -Pi..Pi.
We see that is exactly one (real) solution for any given q in that open interval, not two.

Notice that if q < Pi but close, then the solution x is large.
Try plotting the first found solution res[1] as a function of q:
plot(res[1],q=-Pi..Pi,color=red,discont=true,thickness=3); p1:=%; #Saving the plot in p1
#Observe that the range is rather restricted (in fact to -sqrt(2)..sqrt(2) ). Therefore res[1] will not give the correct answer if q is close to Pi since x is rather large. The second expression will, however.
Try plotting that:
plot(res[2],q=-Pi..Pi,-10..10,color=blue,discont=true,thickness=3); p2:=%: #Saving this plot in p2
##Now show them together:

## The blue-red-blue curve going through (0,0) is the relevant part!
##Confusing? Yes!
### Finally, the good news: For concrete values of q, solve finds the correct formula of the two:
solve(eval(eq,q=3),x);
eval(res[2],q=3);
evalf(%);
solve(eval(eq,q=1),x);
eval(res[1],q=3);
evalf(%);

@Christopher2222 What is shown in the image you bring from mit.edu couldn't be using the OP's U, nor could it be your V:

U:=-0.999047/sqrt((x-(-0.000953))^2 +y^2 )-0.000953/sqrt((x-0.999047)^2 +y^2 );
V:=-10/sqrt((x+10/11)^2 +y^2 )-1/sqrt((x-100/11)^2 +y^2 );

That image has two closed contours not enclosing the singularities. Thus there must inside each of these be either a local minimum or a local maximum. None such exists for U nor V:

solve({diff(U,x)=0,diff(U,y)=0},{x,y});
  Output:   {x = .9690869114, y = 0.}
which is actually the saddle expected to be between the 2 singularities on the y-axis.
The same for your V:

solve({diff(V,x)=0,diff(V,y)=0},{x,y},explicit);
evalf(%);
   {x = 6.688378356, y = 0.} again the saddle point to be expected.

You may want to verify that solve didn't miss any solutions, but that is easily done since the equation diff(U,y)=0 or diff(V,y)=0 clearly only is zero if y=0. Then look at diff(U,x)=0 or diff(V,x)=0 evaluated at y=0 and examine (e.g. plot).

@NomNomPancake I just tested the code in Maple 2015.2. It works.

@NomNomPancake I just tested the code. I had no problems with it in Maple 2015.2.

After having changed the insipid "the season" to "Christmas", I liked it!

@Johnluo In this case at least a numerical ode solution is readily available for comparison:

restart;
myeq:=f(x)=1-int((m(x)*5/3+n(y))*f(y),y=0..x);
ode1:=diff(F(x),x)=f(x);
subs(int(f(y),y=0..x)=F(x),IntegrationTools:-Expand(myeq));
ode2:=subs(ode1,diff(%,x));
m:=x->-0.028+0.0318*(1-exp(-x/10))/(x/10)+0.0657*((1-exp(-x/10))/(x/10)-exp(-x/10));
n:=y->0.0682-0.02*(1-exp(-y/2.479))/(y/2.479)-0.0606*((1-exp(-y/2.479))/(y/2.479)-exp(-y/2.479));
res:=dsolve({ode1,ode2,f(0)=1,F(0)=0},numeric);
plots:-odeplot(res,[x,f(x)],0..15,color=black,thickness=2); p0:=%: #Saving the plot in p0
plots:-odeplot(res,[x,f(x)],0..100,color=black);
##Comparing with:
intsolve(myeq,f(x),method=Neumann,order=1);
res1:=value(%) assuming x>0;
res2:=intsolve(myeq,f(x),method=Neumann,order=2);
plot(rhs~([res1,res2]),x=0..15); p1:=%:
plots:-display(p0,p1);
#From below and up it is order 1, order 2, and the ode solution:

@Johnluo I don't think that you will have any chance with using intsolve for this.

Just consider method=Neumann with orders 1, 2, and 3:

restart;
m:=x->-0.028+0.0318*(1-exp(-x/10))/(x/10)+0.0657*((1-exp(-x/10))/(x/10)-exp(-x/10));
n:=y->0.0682-0.02*(1-exp(-y/2.479))/(y/2.479)-0.0606*((1-exp(-y/2.479))/(y/2.479)-exp(-y/2.479));
myeq:=f(x)=1-1*int((m(x)/0.6+n(y))*f(y),y=0..x);
intsolve(myeq,f(x),method=Neumann,order=1);
res1:=value(%) assuming x>0; #Explicit
res2:=intsolve(myeq,f(x),method=Neumann,order=2); #Explicit
plot(rhs~([res1,res2]),x=0..15);
res3:=intsolve(myeq,f(x),method=Neumann,order=3);
value(%) assuming x>0;
##Output:   f(x) = (x-Float(infinity))/x
## If it could be trusted (surely no) that would be the end. But even if it is wrong there is no hope.
In particular since your actual functions m and n might be much worse than what you gave us.

Notice that intsolve doesn't solve numerically. You need a numerical solver.

There was a similar problem in the link
http://www.mapleprimes.com/questions/204788-Recursive-Integral-Equation

Surely the approach I gave there will work here too.

@Rouben Rostamian  Yes, and that the solution is unique can be seen easily by turning the problem into an initial value problem for a system of odes.

We don't need the special features of m and n besides the fact that they are both continuous (in fact arbitrarily often differentiable), so I don't assign to n or m.

restart;
myeq:=u(x)-int((m(x)+n(y))*u(y),y=0..x);
ode1:=diff(U(x),x)=u(x);
subs(int(u(y),y=0..x)=U(x),IntegrationTools:-Expand(myeq))=0;
ode2:=subs(ode1,diff(%,x));
##This linear ode problem {ode1,ode2,u(0)=0,U(0)=0} has a unique solution by standard uniqueness theorems, and it is obviously as pointed out u(x)=0 and U(x)=0.
#dsolve agrees:
dsolve({ode1,ode2,u(0)=0,U(0)=0},{u(x),U(x)});

@Dayana I ended my reply about the existence of 2 solutions by saying:
"If this is a problem from some applied area, then the second solution is probably the desired one. I'm guessing that N=20 is a replacement for N=infinity."

So is the interval eta=0..20 a replacement for eta=0..infinity?
If so (as I also said), the solutions that becomes (roughly) constant for large eta are most likely the acceptable ones.
If that is the case then it is a waste of time to use as large an interval as eta=0..20 as noted by I_Mariusz.

If it is just a mathematical interest in the possibility of the existence of multiple solutions of this boundary value problem then that is already shown for n=1. By the following code I was able to go as far as n=1.5 for both solutions.
They are done in tandem, so if one breaks one the next is not computed. Actually n=1.6 goes fine for res1, but doesn't for res2. So the fun stops there.

The code could surely be refined, but that is most likely a waste of time.

restart;
Digits := 15;
with(plots):

mu(eta):=(diff(U(eta),eta)^(2)+diff(V(eta),eta)^(2))^((n-1)/(2)):
Eqn1 := 2*U(eta)+(1-n)*eta*(diff(U(eta), eta))/(n+1)+diff(W(eta), eta) = 0;
Eqn2 := U(eta)^2-(V(eta)+1)^2+(W(eta)+(1-n)*eta*U(eta)/(n+1))*(diff(U(eta), eta))-mu(eta)*(diff(U(eta), eta, eta))-(diff(U(eta), eta))*(diff(mu(eta), eta)) = 0;
Eqn3 := 2*U(eta)*(V(eta)+1)+(W(eta)+(1-n)*eta*U(eta)/(n+1))*(diff(V(eta), eta))-mu(eta)*(diff(V(eta), eta, eta))-(diff(V(eta), eta))*(diff(mu(eta), eta)) = 0;
bcs1 := U(0) = 0, V(0) = 0, W(0) = 0;
bcs2 := U(N) = 0, V(N) = -1;

AP1:=NULL;
AP2:=approxsoln=[U(eta)=eta/20*exp(-eta),V(eta)=-tanh(eta/2),W(eta)=-0.9*tanh(eta/2)];
for n from 1 to 1.9 by 0.1 do
  res1[n]:=dsolve(eval({Eqn1, Eqn2, Eqn3,bcs1, bcs2},N=20), initmesh=10000,numeric,AP1,abserr=1e-5);
  AP1:=approxsoln=res1[n];
  try
    res2[n]:=dsolve(eval({Eqn1, Eqn2, Eqn3,bcs1, bcs2},N=20), initmesh=10000,numeric,AP2,abserr=1e-5);
  catch:
    AP2:=approxsoln=res2[n-0.1];
    res2[n]:=dsolve(eval({Eqn1, Eqn2, Eqn3,bcs1, bcs2},N=20), initmesh=10000,numeric,AP2,abserr=1e-5);
    AP2:=res2[n]
  end try;
  printf("success for n = %a\n",n)
end do;
 
indices(res1);
indices(res2);
odeplot(res1[1.6],[eta,U(eta)]);
odeplot(res2[1.5],[eta,U(eta)]);

@Maple4evergr8 You can obtain separate files (plot1.jpg,plot2.jpg, ...) by doing:

FileTools:-MakeDirectory("E:/MapleTemp");
mypath := "E:/MapleTemp/plot":
for i from 1 to 10 do
    p := plot(x^i,x=-1..1);
    plottools:-exportplot(cat(mypath,i,".jpg"),p)
end do:

There are quite a few tools in the package StringTools:

?StringTools

Continuing where I left above saying that the difference in results merits looking into.

There are at least two solutions for n=1 and N=20, and I think for other values as well.

The nice thing about n=1 is that the system is so simple.

Again for completeness I bring the whole code. I don't assign to n or N, but use eval.

restart;
Digits := 15;
with(plots):
mu(eta):=(diff(U(eta),eta)^(2)+diff(V(eta),eta)^(2))^((n-1)/(2)):
Eqn1 := 2*U(eta)+(1-n)*eta*(diff(U(eta), eta))/(n+1)+diff(W(eta), eta) = 0;
Eqn2 := U(eta)^2-(V(eta)+1)^2+(W(eta)+(1-n)*eta*U(eta)/(n+1))*(diff(U(eta), eta))-mu(eta)*(diff(U(eta), eta, eta))-(diff(U(eta), eta))*(diff(mu(eta), eta)) = 0;
Eqn3 := 2*U(eta)*(V(eta)+1)+(W(eta)+(1-n)*eta*U(eta)/(n+1))*(diff(V(eta), eta))-mu(eta)*(diff(V(eta), eta, eta))-(diff(V(eta), eta))*(diff(mu(eta), eta)) = 0;
bcs1 := U(0) = 0, V(0) = 0, W(0) = 0;
bcs2 := U(N) = 0, V(N) = -1;
sys:=eval({Eqn1, Eqn2, Eqn3},n=1); #Relatively simple
res1:=dsolve(sys union eval({bcs1, bcs2},N=20), numeric); #Using N=20
odeplot(res1,[eta,U(eta)],thickness=3); p1:=%:
odeplot(res1,[eta,V(eta)],thickness=3);
odeplot(res1,[eta,W(eta)],thickness=3);
##Now using an appropriate approximate solution to find the second solution:
res2:=dsolve(sys union eval({bcs1, bcs2},N=20), numeric,
    approxsoln=[U(eta)=eta/20*exp(-eta),V(eta)=-tanh(eta/2),W(eta)=-0.9*tanh(eta/2)]);
odeplot(res2,[eta,U(eta)],thickness=3,color=blue); p2:=%:
odeplot(res2,[eta,V(eta)],thickness=3,color=blue);
odeplot(res2,[eta,W(eta)],thickness=3,color=blue);
display(p1,p2,caption="Two solutions for U");


If this is a problem from some applied area, then the second solution is probably the desired one. I'm guessing that N=20 is a replacement for N=infinity.

Comment. If you use bcs3:=D(U)(N) = 0, D(V)(N) = 0;  instead of bcs2 then the red solution is excluded.

@Dayana It seems that the statement in the help page about initmesh is incorrect:

?dsolve,numeric,bvp

saying (quoting the entire statement):

" 'initmesh'= integer
Integer value that determines the number of points dsolve uses to compute the initial solution profile. In some cases, the default initial 8 point mesh does not have sufficient resolution to obtain the initial solution profile, so increasing this value can give a solution when the default value does not. Its value must be between 8 and 8192."
(My emphasis)

Consider the following simple bvp-problem (whose solution is just y(x) = cos(x) ).
Compare the different usages reported.

restart;
ode:=diff(y(x),x,x)+y(x)=0;
res:=CodeTools:-Usage(dsolve({ode,y(0)=1,D(y)(Pi)=0},numeric,initmesh=30000));
     memory used=86.52MiB, alloc change=83.83MiB, cpu time=3.67s, real time=3.67s, gc time=46.88ms
restart;
ode:=diff(y(x),x,x)+y(x)=0;
res:=CodeTools:-Usage(dsolve({ode,y(0)=1,D(y)(Pi)=0},numeric,initmesh=8192));
     memory used=26.82MiB, alloc change=47.97MiB, cpu time=516.00ms, real time=516.00ms, gc time=0ns
restart;
ode:=diff(y(x),x,x)+y(x)=0;
res:=CodeTools:-Usage(dsolve({ode,y(0)=1,D(y)(Pi)=0},numeric)); #Default
     memory used=4.87MiB, alloc change=32.00MiB, cpu time=47.00ms, real time=46.00ms, gc time=0ns

I think that we must conclude that the statement "Its value must be between 8 and 8192." is incorrect at least all the way back to Maple 12, which is the oldest I have available at my present location.

Also I shall submit an SCR asking for a help page correction (update).
There is a similar statement about maxmesh, but doing the same experiment with maxmesh instead of initmesh doesn't show any significant difference (try several times). That is to be expected since maxmesh is the maximal mesh allowed and if the mesh doesn't need to be raised above the given maxmesh then it isn't.

Thanks for your persistence in your inquiries!

@Dayana When getting the error 'initial Newton iteration not converging' you could try using the option 'initmesh'= N, where N is between 8 and 8192. (8 is the default according to the help for dsolve,numeric,bvp).

Another useful option is to give an approximate solution in the form 'approxsoln'= [U(x)=..., V(x)= ..., W(x)=...].
You can try any expression depending on (in this case) eta.

A third option to consider is the absolute error in the form 'abserr'= eps, where eps is some small positive number. The default is 10^(-6).

The continuation method as you mention is more difficult to use because it involves placing a parameter in a "reasonable" place.

As far as using a midpoint method is concerned, dsolve will often complain itself about problems at an endpoint and suggest using a midpoint method.