@shahid Since you have a Maple version earlier than Maple 18 eval[recurse] doesn't work.

This should work (I did this in Maple 15):

restart;
eq1:=u(x,y)*(diff(u(x,y), x))+v(x,y)*(diff(u(x,y), y)) = (diff(u(x,y), y, y))/(1+epsilon*theta(x,y))-epsilon*(diff(u(x,y), y))*(diff(theta(x,y), y))/(1+epsilon*theta(x,y))^2+theta(x,y);

eq2:=u(x,y)*(diff(theta(x,y), x))+v(x,y)*(diff(theta(x,y), y)) = (diff(theta(x,y), y, y))/Pr;

## Replace the original S with the following:
S1:={u(x,y) = diff(psi(x,y), y), v(x,y) = -(diff(psi(x,y), x))};
S2:={ psi(x,y) = x^(4/5)*f(x,eta)/(1+x)^(1/20), theta(x,y)=x^(1/5)*(1+x)^(-1/5)*h(x,y)};
S3:= eta = y/(x^(1/5)*(1+x)^(1/20));
S23:=eval(S2,S3);
eval([eq1,eq2],S1);
eval(%,S23);
subs(rhs(S3)=lhs(S3),%);
collect((lhs-rhs)~(%)=~0,[D,diff],factor);
res:=subs(y/(x^(1/5)*(1+x)^(1/20))=eta,%);

##To make it more readable you can continue with:
PDEtools:-declare(f(x,eta),h(x,y));
convert(res,diff);

About rewriting S1*S2-S3*S4.

You can rewrite the products of the inert sums S1*S2 and S3*S4 using the Cauchy product. Thus you get

S1*S2-S3*S4 = S12 - S34, where S12 is the Cauchy product of S1 and S2, and S34 is the Cauchy product of S3 and S4.

The difference S12-S34 can be combined to just one sum S:

S := 1/(2*Pi)*Sum(x^(3*k1+2)*3^(1/2+k1)*(Sum((3*GAMMA(k+2/3)*(-2/3-k1+k)*GAMMA(k1+1/3-k)+3*GAMMA(4/3+k)*GAMMA(2/3+k1-k))/(GAMMA(3*k+2)*GAMMA(3*k1-3*k+3)), k = 0 .. k1)), k1 = 0 .. infinity);

This can be written in short form as
1/(2*Pi)*Sum(x^(3*k1+2)*3^(1/2+k1)*m(k1), k1 = 0 .. infinity);

where m(k1) < 0 for all k1 >=0. It is given by:

m(k1) = Sum((3*GAMMA(k+2/3)*(-2/3-k1+k)*GAMMA(k1+1/3-k)+3*GAMMA(4/3+k)*GAMMA(2/3+k1-k))/(GAMMA(3*k+2)*GAMMA(3*k1-3*k+3)), k = 0 .. k1);

m(k1) is very rapidly decreasing in absolute value.

Now you don't have cancellation problems for x > 0 since all terms in S are negative.
But I wonder why x-values as high as x=100 are interesting to you. S will be huge (but negative) for x = 100.
For x=10 you get
evalf[50](eval(S,x=10));
    -4.7714630011989009659113681348806733279336141633480*10^8
and for x = 11:
evalf[50](eval(S,x=11));
    -1.1891747657696851151848674770800968156039786693327*10^10
### Notice that you will run into similar problems with the well known power series for exp(x):
restart;
S:=Sum(x^n/n!,n=0..infinity);
seq(evalf[500](eval(S,x=i*10)),i=0..10);

@shahid Here is what I got (last output shown only) when I just now copied my own code from above and executed the whole thing in Maple 2015.2. Actually, in my worksheet it is shown with the abbreviated notation used because of my use of PDEtools:-declare(f(x,eta),h(x,y)). But when I copy and paste it into MaplePrimes it shows up as you see below.

You could begin by telling us what Ni1 is. Maybe you also can tell us if a1 and a2 are assigned to anything.

@Bendesarts You wrote:

"Don't hesitate to let me know if you have new ideas for troubleshooting my issue of "drift" for the solution of this NL oscillator but in MapleSim."
(my emphasis)

I don't have MapleSim.

@shahid You forgot theta(x,y)=x^(1/5)*(1+x)^(-1/5)*h(x,y) it seems. I added it to S below.

restart;
eq1:=u(x,y)*(diff(u(x,y), x))+v(x,y)*(diff(u(x,y), y)) = (diff(u(x,y), y, y))/(1+epsilon*theta(x,y))-epsilon*(diff(u(x,y), y))*(diff(theta(x,y), y))/(1+epsilon*theta(x,y))^2+theta(x,y);

eq2:=u(x,y)*(diff(theta(x,y), x))+v(x,y)*(diff(theta(x,y), y)) = (diff(theta(x,y), y, y))/Pr;

S:={ psi(x,y) = x^(4/5)*f(x,eta)/(1+x)^(1/20), eta = y/(x^(1/5)*(1+x)^(1/20)), u(x,y) = diff(psi(x,y), y), v(x,y) = -(diff(psi(x,y), x)),theta(x,y)=x^(1/5)*(1+x)^(-1/5)*h(x,y)};

eval[recurse]([eq1,eq2],S);
collect((lhs-rhs)~(%)=~0,[D,diff],factor);
res:=subs(y/(x^(1/5)*(1+x)^(1/20))=eta,%);

##To make it more readable you can continue with:
PDEtools:-declare(f(x,eta),h(x,y));
convert(res,diff);

@Bendesarts Did you try to run your code in Maple itself, not only in MapleSim?

Clearly, we don't know (or rather I don't) whether there is actually a limit cycle. But at least it appears (or appeared) so.
You have an autonomous system of 8 odes of first order. If a given orbit stays bounded, the limit set may be quite complicated, just think of the Lorenz system of only 3 odes.
Even for a system of just 2 odes (where the possibilities for the limit set are well known) it may be difficult to determine numerically whether you have a limit cycle or just a very, very slow oscillating approach to an equilibrium point.

You may want to change the name of that attached file.

@shahid You should present your equations or expressions using Maple code as you did in your original question.

@Athar Shahabinejad You ought to tell us what it is you are running now, when you say

"Now when i run i get this answer ... "

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