@Christian Wolinski To answer your question:

Why is the output not:
['a', 1 .. 4, x, 'b', 1 .. 4, x, 'c', 1 .. 4, x, 'd', 1 .. 4, x]
and
[a, '1 .. 4', x, b, '1 .. 4', x, c, '1 .. 4', x, d, '1 .. 4', x]

my reply is that you should not expect quotes in the output at all.
Consider the more elaborate version of the first procedure (a::uneval ,b)->'args', now using the syntax in Maple V Release 4.
p:=proc(a::uneval,b) 'args' end;
p(x,y,z); # returns x,y,z in all Maple versions, including Maple V Release 4.
This is as it should be since first args is simply replaced by x,y,z and then evaluation to one level removes the quotes around 'x,y,z'. (Full evaluation would do just the same).

Now to illustrate the effect of the declaration a::uneval, let
xx:=89;
# and try
p(xx,y,z); #Output xx,y,z since evaluation of 'xx,y,z' one level just means removing the unevaluation quotes. (Again: Full evaluation would do just the same, but we are inside a procedure).
# If you continue with
"; # That is % these days.
#  then you get 89,y,z
Now consider instead
q:=proc(a::uneval,b) 'args[1..-1]' end;
q(xx,y,z); #Output (xx,y,z)[1..-1] as it should.
# Again if you continue with
"; # You get full evaluation to 89,y,z

@necron I just tried Digits=100. That helped, but didn't produce the graph you expected.

To your PS: It is even harder to help if there is no worksheet. You gave us images, thus we would have to type the whole thing ourselves.

@torabi With your new data as presented in your worksheet DSOLVE3.mw I had success with using continuation as follows.

1. Proceed from restart until you have defined sys.
2. Then do the following:

indets(sys,specfunc(ln)); #Just to see the culprits
##We try making the ln terms harmless:
sys0:=evalindets(sys,specfunc(ln),x->ln(10));
res0:=dsolve({op(sys0)} union bcs2 union {bcs22}, numeric);
plots:-display(Array([seq(plots:-odeplot(res0,[theta,h(theta)]),h=[h1,h2,h3])]));
## That was fine. Based on that experience we introduce a continuation parameter lambda
## such that lambda=0 gives us sys0 and lambda=1 gives us sys:
sysC:=evalindets(sys,specfunc(ln),x->ln(subs(h3(theta)=lambda*h3(theta),op(1,x))));
## Now use continuation in lambda:
res:=dsolve({op(sysC)} union bcs2 union {bcs22}, numeric,continuation=lambda);
plots:-display(Array([seq(plots:-odeplot(res,[theta,h(theta)]),h=[h1,h2,h3])]));

@torabi In what way does the new version differ from the old. You are using new data, you say. How far are they from the data in the first version? Could you use the solution from that for an approximate solution?
The error message probably is due to the log terms becoming imaginary.

I shall be gone for a week so won't have time for more now.

@tomleslie In his original version, which I used in my answer, I get the output 3 from both of

nops(SYS);
nops(SYS2);

This reminds me of a similar question:

http://mapleprimes.com/questions/215244-Unevaluated-Expression#answer230511

@Carl Love  I saw some examples by Axel Vogt about this, but I don't have the reference.
Here is a contrived example made by me, but I think it catches one of his main points:
kernelopts(floatPi=false);
seq(sin(2*n*Pi+0.5),n=0..10000,1000); #Exactly the same
kernelopts(floatPi=true);
seq(sin(2*n*Pi+0.5),n=0..10000,1000); #Varies as can be expected

## Of course you could make this example more glaring if you began with e.g.  Digits:=5:

@Carl Love I think my primary reason is consistency with the treatment of other constants:

kernelopts(floatPi=true);
Pi+2.;         #float
Catalan+2.; #mix
gamma+2.;  #mix
exp(1)+2.;   #mix
sin(2)+2.;    #mix

@Carl Love It may be relevant to mention that since Maple 2015 the output of input like

2.0 + Pi;

can be controlled by setting
kernelopts(floatPi=true);
or
kernelopts(floatPi=false);

The option is available (but not mentioned) in Maple 2015. Edited:  The default is true, i.e. output is 5.141592654, whereas it is 2.0+Pi in all previous versions of Maple.
The option is available (and mentioned) in Maple 2016. The default is true, i.e. output is 5.141592654.

I have put kernelopts(floatPi=false);  in my maple.ini file.
Added: Actually, because kernelopts(floatPi=false) will generate an error in Maple 18 and earlier versions i have wrapped the statement in  try catch end try:

try kernelopts(floatPi=false) catch: end try:

@nm Actually the help page on pi (?pi) has as its first line:

"For information on 3.14159..., see Pi."

Furthermore, the two relevant examples given on that page uses Pi.

Although I never use 2D math input myself I have noticed that in that context there is some confusion (second-guessing).
Example:
In a fresh worksheet using 1D input (Maple input) I write

a:=pi

Then I select that and right click to bring up a context menu.
I choose Convert To/ 2-D Math Input.
The response is  

as expected because pi and Pi print the same.
When I use
lprint(a);
I get Pi, and when I use
evalf(a);
I get  3.141592654.

@mehdi jafari So Xstar is a function of T. Thus the first equation should be written maybe as
EQ[1]:=(1/Dg)*diff(X(T,Z),T)+diff(X(T,Z),Z)+3*St*(X(T,Z)-Y(T,1)^(1/n));

I don't see any way to handle this equation together with EQ[2].

@Carl Love I agree. As I said in my answer: If you start at (theta,theta_p) =(Pi,0) you will never leave, if you don't, you will never get there!
So Kitonum's remark in his answer: "At this value there is instability when the pendulum reaches the top position." is misleading. The pendulum cannot reach the top position and also be at rest!

@mskalsi I agree that frontend is not easy to understand, so I appreciate the freeze/thaw approach used by Carl.

Incidentally, the first term in your input expression had the factor diff(g(z),x), where surely you meant diff(g(z),z).

To give an example, where more freezing is necessary (and by default handled by frontend) take this modified expression:

eq2 := sin(g(z))*diff(g(z), z)*g(z)^3+diff(g(z), z, z)*g(z)^4+diff(g(z), z, z, z)*g(z)^5+diff(g(z), z)/g(z)^2;

The command
thaw(coeff(subsindets[flat](eq2, specfunc(diff), freeze), g(z),3));
#doesn't succeed because sin(g(z)) won't be frozen and because coeff doesn't work here:
coeff(sin(g(z))*F*g(z)^3,g(z),3);
#By default frontend also freezes sin(g(z)) as is illustrated here
frontend(diff,[sin(g(z)),z]); #returns zero
frontend(diff,[sin(g(z)),z],[{`+`,`*`,function},{}]); #Doesn't freeze expressions of type function.
#So frontend succeeds in
frontend(coeff,[eq2,g(z),3]);
# but not in
frontend(coeff,[eq2,g(z),3],[{`+`,`*`,function},{}]);

@Kitonum It appears that it is the cook method which gives erroneous results:
restart;
res:=[seq(int(1/(x^4+p), x=0..1,method=_RETURNVERBOSE), p=1/2..5, 1/2)];
evalindets(res,`=`(identical(FAILS),anything),x->NULL);
evalf(%);
restart;
res:=[seq(int(1/(x^4+p), x=0..1,method=cook), p=1/2..5, 1/2)];
###################################
I tried for general p with and without the assumption p>0.

## Without assumptions:
restart;
infolevel[int]:=3:
res:=int(1/(x^4+p), x=0..1);
simplify(evalc(eval(res,p=4))); #cook correct
evalf(%);
##With the assumption p>0:
restart;
infolevel[int]:=3:
res:=int(1/(x^4+p), x=0..1) assuming p>0;
simplify(eval(res,p=4)); #cook wrong

## Just in case this bug has not been reported as an SCR, I shall do it.