@acer Yes, the URL for that discussion is
http://www.mapleprimes.com/posts/201844-Possible-Bug-In-Maple-

The problem seems to be that I had the following maple.ini file, which solved some other problem. This fix was supplied by Edgardo ChebTerrab (and acer), but that I used it I cannot blame on them.

__fixplus:=module() option package; export `+`; end module:

`print/+`:=eval(:-`+`):

with(__fixplus):

macro(__fixplus:-`+`=:-`+`):

The Lobatto procedure supplied by Carl works as written without this redefinition of `+`.
All my examples above work.

@Carl Love Please see correction below.

You are quite right. In my opinion this is a bug.
It seems that D doesn't like expressions containing indeterminates of type `+`  (!)

I have often and for many years been disappointed with D as I had thought that it would at least be able to do the work that diff could, but I quickly learned that that is far from the case.
This may be due to a misunderstanding on my side of the purpose of D, but if so, the limitations ought to be stated clearly in the help page. The help page contains an example:

f := (x,y) -> exp(x*y);
D[1](f);
# That example would break, by just adding a 1:
f := (x,y) -> exp(x*y)+1;


More examples:

restart;
f:=x->x^2*sin(x);
D(f); #OK
D(f)(Pi/2); #OK
(D@@2)(f); #OK
(D@@2)(f)(Pi/2); #OK
##
f:=x->x^2+sin(x); #A sum
indets(f(x),`+`);
D(f); #No evaluation
D(f)(Pi/2); #No evaluation
evalf(%); #Numerical result OK
evalf((D@@2)(f)(Pi/2)); #No evaluation
##
f:=x->x^2-1;
D(f); #No evaluation
D(f)(2); #No evaluation
evalf(%); #OK
evalf((D@@2)(f)(2)); ##No evaluation
##
f:=x->2^(x+1);
indets(f(x),`+`);
D(f); #No evaluation
D(f)(-1); #No evaluation
evalf(%); #OK
evalf((D@@2)(f)(-1)); ##No evaluation

I shall submit an SCR although it must be well known.

## An additional weird thing:
restart;
f:=x->x^2*sin(x); #My first example
D(f); #OK
(D@@2)(f); #OK
D(f); #OK. The result is a sum, so maybe the next ought not work:
D(%); #But it does! But now try:
f1:=x->2*x*sin(x)+x^2*cos(x);
D(f1); #No evaluation
D(x->2*x*sin(x)+x^2*cos(x)); #No evaluation

@Carl Love I tried your Lobatto procedure in your worksheet on the example you have. I got



NOTE: This result was due to a redefinition of `+` I had in a maple.ini file. That redefinition was supposed to solve some other problem. That redefinition apparently spoils D.

Changing the procedure to the following helped. I wonder why it worked for you? I used Maple 2015.2.

Lobatto:= proc(f::algebraic, R::name=range(realcons), n::{posint, Not(identical(1))})
local x, a, b, P,DP,F,oldDigits,r;
     x:= lhs(R);
     (a,b):= op(evalf(op(2, R)));
     P:= unapply(diff((x^2-1)^(n-1)/2^(n-1)/(n-1)!,x$(n-1)),x); ##The important part
     DP:= D(P)(x);
     oldDigits:= Digits;
     Digits:= Digits+1+ilog2(Digits)+iroot(n,3)^2;
     F:= unapply(eval((b-a)*f, x= (b-a)/2*x + (a+b)/2)/P(x)^2/n/(n-1), x);
     r:= (b-a)*(eval(f, x= a)+eval(f, x= b))/n/(n-1) + add(F(x), x= [fsolve(DP)]);
     evalf[oldDigits](r)
end proc:

@shadi1386 Do it numerically:

restart;
eq:=(19.42795980*(1+z))/x(z)-(14.09324088*(1+z))*(eval(diff(x(z), z), z = 0))/(x(z)*x(0))-19.42795980*(1.+z)^(3/2)+14.09324088*(diff(x(z), z))/(x(z)*(1/(1+z)^(3/2))^(5/3))-19.42795980-19.42795980*z+19.42795980*sqrt(.314*(1+z)^3+.686)=0;
ode:=subs(x(0)=x0,eval(diff(x(z), z), z = 0)=x1,eq);
eval(ode,z=0);
eq0:=eval(%,{x(0)=x0,eval(diff(x(z), z), z = 0)=x1});
x11:=solve(eq0,{x1});
ode1:=eval(ode,x11);
res:=dsolve({ode1,x(0)=2},numeric);
plots:-odeplot(res,[z,x(z)],-0.8..3);
res:=dsolve({ode1,x(0)=1},numeric);
plots:-odeplot(res,[z,x(z)],-0.8..3);

@Mac Dude I get from the exact lines I had in Maple 2015.2:



But I realize that I deleted your use of the Physics package as it looked completely horrible (no exaggeration) in the 1D notation you supplied in your worksheet.

##Note: I just tried this:

restart;
with(Physics[Vectors]):
Setup(mathematicalnotation=true);
E:=-2*P1*A[Q1]+A[0]^2+A[Q1]^2-2*A[0]*rho*P2/(rho+Q1)+A[Q3]^2-2*P3*A[Q3]+P1^2+rho^2*P2^2/(rho+Q1)^2+P3^2+m^2*c^4;
mtaylor(E,[A[Q1]=P1,A[0]=P2*rho/(rho+Q1),A[Q3]=P3],6);
Student:-Precalculus:-CompleteSquare(E,[A[0],A[Q1],A[Q3]]);

and the result is the same. So I don't know what is going on at your end,

@Markiyan Hirnyk When you "'manipulate" an equality or inequality by using manipulators as e.g. map, subsop, evalindets, or subs, you are the one responsible for the correctness. The manipulators will just do the manipulation you asked for.

Trivial example:
eq:=a+b=c;
evalindets(eq,name,sin);

@Bendesarts What I see in your MapleSim plot is a curve approaching the limit cycle, which seems to be there.
Try repeating the run with the the same initial values at 0, but showing only the curve for t=8..tmax.
I don't have MapleSim, so I cannot try myself, but in Maple I did try the 3d plot:
plots:-odeplot(res,[u[1](t),v[1](t),u[2](t)],8..100, refine=1,scaling = constrained);

where I (rather arbitrarily) used u[2](t) as the third variable.


Incidentally, you can optionally use method=ck45 in dsolve, if you wish. Also relerr and abserr can be given as 1e-5 as in your MapleSim image. That doesn't change the appearance of the picture I show above.

@Bendesarts The movement should (very nearly) be periodic with the period of the limit cycle. The period of the limit cycle could be determined like this (this for your last example):

res:=dsolve([sys[],ic[]],numeric, range=0..8,output=listprocedure);
U1:=subs(res,u[1](t));
u1_8:=U1(8);
plot(U1,7..8);
fsolve(U1-u1_8,[8-0.6]);
T:=8-%;
# I get 0.558282997

I don't see any point in using dsolve all the way to t=300 since you can just use the periodicity obtained rather early (here before t=8). Of course there is roundoff, but what about your actual robot? Is that ideal?

Maybe I don't understand your problem?

@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!