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

@torabi  Just a comment: Your latest system is linear, so there is no particularly good reason to use small values of b like b=10^(-5). You might as well use b=1.

@ThU Indeed, Maple 17 works on my 64 bit computer with Windows 10 Home Edition.

@Carl Love To put the system into a form that allows you to solve for the highest derivatives you can differentiate the first of the equations. This won't introduce any higher derivatives in the system:

ode1:=diff(sys[1],x);
solve({sys[2],ode1}, {diff(w(x),x$4), diff(psi(x),x$3)}); # Possible
# Now add a boundary condition to ensure that the new system is equivalent to the old:
bcs1:=eval(convert(sys[1],D),x=1); # Might as well have used x = 0.
# So new boundary conditions are
BCS:=bcs union {bcs1};
nops(BCS);
## To these must be added the extra condition to determine omega2.
##
By doing that I get 8 successes out of 9 tries without any optional arguments to dsolve expect for the output option, i.e. only
dsolve({sys[2],ode1} union BCS union {b = 10^(-5)}, numeric,output=listprocedure) . 

'undefined' is not the result of an error.
Example:
Heaviside(0);
                                  undefined
But 0/0 produces the error
                            Error, numeric exception: division by zero
i.e. it doesn't result in 'undefined'.

I suppose that you have already looked at the help page:

?undefined

@mwahab I don't know much about this (actually close to nothing), but I tried the following where my result RES differs from the result obtained by map(pdsolve, [res2], parameters = {m}) in the latter missing the case m = 0.

restart;
with(PDEtools);
declare(u(t, x))
pde:=diff(u(t, x), t)-u(t,x)^m*(diff(u(t, x), x))-u(t,x)^m-u(t,x)^m*(diff(u(t, x), x, x))-u(t,x)^m*(diff(u(t, x), x, x, x)) = 0;
res:=DeterminingPDE(pde, u(t, x), integrabilityconditions = false);
nops(res);
res2:=casesplit(res, parameters = {m});
map(pdsolve, [res2], parameters = {m});
n:=nops([res2]);
map2(op,2,[res2]); # The &where's
eqs:=map2(op,1,[res2]); #The equations including conditions on m
sys:=map2(select,has,eqs,{_eta,_xi}); #the pdes
M:=map2(remove,has,eqs,{_eta,_xi}); # the m equations
M1:=subs([]=[m=m,m=m],M); #Slightly artificial
M2:=map2(op,2,M1); #Final version of m equations
RES:=seq( [M2[i],pdsolve(eval(sys[i],M2[i]))],i=1..n); #result
map(pdsolve, [res2], parameters = {m});

See the help page for dsolve:
?dsolve

Include any initial conditions to determine the arbitrary constant(s).

Have a look at the help page
?dsolve,numeric

It has examples of ivp as well as bvp problems.

Note: I must agree with vv that this seems to be a case of using a generic formula.
So I converted this into a comment.
Here is a simplified version, where I set t = 0 and a=b=2 right away.

restart;
G:=n->sum(cos(2*Pi/n*j)^2*sin(2*Pi/n*j)^2,j=0..n-1); #The simplified version
Ga:=n->add(cos(2*Pi/n*j)^2*sin(2*Pi/n*j)^2,j=0..n-1); #Version using add
H:=n->sum(combine(cos(2*Pi/n*j)^2*sin(2*Pi/n*j)^2),j=0..n-1); #Version combining before summing
G(n); #Wrong for some n
H(n); #Seems to be better
combine~([seq(G(n),n=1..5)]);
combine~([seq(Ga(n),n=1..5)]);
[seq(H(n),n=1..5)];

##Test of n/8 as the answer:
simplify(combine~([seq(G(n)-n/8,n=1..25)]));
simplify(combine~([seq(Ga(n)-n/8,n=1..25)]));
simplify([seq(H(n)-n/8,n=1..25)]);