Maybe frontend is what you need:

eq := (diff(g(z), x))*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;
frontend(coeff,[eq,g(z)^4]);
frontend(coeff,[eq,g(z)^(-2)]);

Although equality assumptions are allowed (I think primarily for reasons of consistency), they don't work as evaluation.
Simple example:
restart;
sqrt(a^2) assuming a=2;

The output is a not 2.
Thus the properties the integer 2 has (like being positive) have been used, but a has not been replaced by 2.

Set Digits higher than the default whish is 10, doing like this: Digits:=50: (or whatever you like!).
Alternatively do
restart;
evalf[50](BesselI(0,31.7)); # or use another integer than 50, e.g. 15
## I get
                 4.1615582165614032637527525720570080615945783690515*10^12
# Notice 10^12, not 10^20.
# Similarly
evalf[50](BesselI(0,60));
# gives me  
       5.8940770556098011682788174403339047379789830203162*10^24

First a theoretical remark:
Rewrite your ode as a system:
DEtools[convertsys](ode,{},theta,t,Y,YP);
   [[YP[1] = Y[2], YP[2] = sin(Y[1])*cos(Y[1])-9.8100*sin(Y[1])], [Y[1] = theta(t), Y[2] = diff(theta(t), t)], undefined, []]
You have an autonomous system in Y= (Y[1],Y[2]). The system has an equilibrium at Y = (Pi,0).
If you start at (Pi,0) you will never leave, if you don't, you will never get there!
So the short answer is: No solution!
##Further comments:
The equilibrium (Pi,0) is a saddle point, thus unstable.
The jacobian (where I use g for 9.81):
J:=unapply(VectorCalculus:-Jacobian([Y[2],sin(Y[1])*cos(Y[1])-g*sin(Y[1])],[Y[1],Y[2]]),[Y[1],Y[2]]):
J(Pi,0);
LinearAlgebra:-Eigenvalues(%); #One positive and one negative, so a saddle point.
A saddle point has 4 separatrices corresponding to solutions satisfying Y(t) -> (Pi,0) as t -> infinity or as t->-infinity.
The directions of approach near (0,Pi) are given by the eigenvectors.
Below we find a solution which gets very close to a separatrix. The solution is approaching (0,Pi) and gets very close, but leaves the vicinity again. See the animation.

So try this:

restart;
ode:=sin(theta(t))*cos(theta(t))-9.8100*sin(theta(t))-(diff(theta(t), t, t)) = 0;
sol:=dsolve({ode, theta(0) = 2*Pi*(1/3), D(theta)(0) = th1}, numeric,parameters=[th1],output=listprocedure);
TH,TH1:=op(subs(sol,[theta(t),diff(theta(t),t)]));
## The way parameters are set:
sol(parameters=[1.0341*Pi]); #Example
plots:-odeplot(sol,[[t,theta(t)-Pi],[t,diff(theta(t),t)]],t=0..5) ;
##Now make a procedure with output we want to be [0,0]:
p:=proc(th1,t0) option remember; sol(parameters=[th1]); [TH(t0)-evalf(Pi),TH1(t0)] end proc;
## Because of LSSolve (or fsolve) we need these:
p1:=proc(th1,t0) p(_passed)[1] end proc;
p2:=proc(th1,t0) p(_passed)[2] end proc;
## Testing:
p1(1.0341*Pi,1.5);
p(1.0341*Pi,1.5);
# I haven't had success with fsolve, so I use LSSolve:
res:=Optimization:-LSSolve([p1,p2],initialpoint=evalf([1.0340*Pi,1.5]));
res[2][1]/evalf(Pi);
plots:-odeplot(sol,[[t,theta(t)-Pi],[t,diff(theta(t),t)]],t=0..6) ;
## Phase plane:
plots:-odeplot(sol,[[theta(t),diff(theta(t),t)]],t=0..10,frames=50,caption="Very close to a separatrix") ;

################## Finding an equation for the separatrices:
ode:=sin(theta(t))*cos(theta(t))-9.81*sin(theta(t))-(diff(theta(t), t, t)) = 0;
## ode implies that the following is zero:
lhs(ode)*diff(theta(t),t);
#So we have
eq:=int(%,t)=C; #Where C is a constant
#Finding C corresponding to a separatrix for (Pi,0):
eval[recurse](convert(eq,D),{t=0,theta(0) = Pi, D(theta)(0)=0});
#So an equation for a separatrix is
eq2:=subs(diff(theta(t),t)=theta_p,theta(t)=theta,eval(eq,C=-9.81));
plots:-implicitplot(eq2,theta=0..15,theta_p=-7..7,gridrefine=2,labels=[theta,theta__p]);

As Markiyan is saying the two solutions are not the same for all values of a and b.
With assumptions a>0, b>a, they are:

odetest(y(x)=sol_2,ode);
combine(%) assuming a>0,b>a;
simplify(%) assuming a>0,b>a;
##
dsolve gives the solution on this form:
res:=dsolve(ode);
       
Assuming that a and b are real then this answer is real for a^2>b^2 (taking _C1 real).
If a^2<b^2 then it is still fine but involves complex parts. It can be rewritten, though.
subs(sqrt(a^2-b^2)=I*sqrt(-a^2+b^2),res);
convert(%,tan);
              
The cases b^2=a^2 are special and can be handled separately, here e.g. the case b=-a:
ode1:=diff(y(x),x)=a*cos(y(x))-a;
where the answer is found as
y(x) = 2*arctan(1/(a*(_C1+x)))

That you get a response like `[Length of output exceeds limit of 1000000]` is not an error.
Just try saving the result:

res := dsolve(dsys3, numeric, initmesh = 3024, abserr = 0.1e-4);

Then try e.g.
res(0.1e-6);

Unsurprisingly you get zeros across the board (except for the value of x).

Clearly your bvp problem has the trivial solution u(x) = 0, w(x) = 0 for all x in the interval between the end points.
Do you have reason to believe that a nontrivial solution exists?

@tomleslie mwahab left out the arguments in u^m in his worksheet, so his pde is

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];

He uses declare(u(t,x)) , which will display u(t,x) as u, but it doesn't allow you to use the shortcut u for u(t,x) in input, at least not always, but unfortunately sometimes! This is not a declare problem though.

Examples.
#The present one.
restart;
pde2 := 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; #Arguments not left out
pde := diff(u(t, x), t)-u^m*(diff(u(t, x), x))-u^m-u^m*(diff(u(t, x), x, x))-u^m*(diff(u(t, x), x, x, x)) = 0; #Arguments left out
pdsolve(pde2); #No answer
pdsolve(pde); #An answer involving u^m, obvious nonsense.
## The same with declare:
restart;
with(PDEtools,declare);
declare(u(t,x));
pde2 := 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; #Arguments not left out
pde := diff(u(t, x), t)-u^m*(diff(u(t, x), x))-u^m-u^m*(diff(u(t, x), x, x))-u^m*(diff(u(t, x), x, x, x)) = 0;#Arguments left out
pdsolve(pde2); #No answer
pdsolve(pde); #Same nonsense as above (but in compact form)
#### A simpler example where leaving out arguments in input works as hoped for (??).
restart;
pde2:=diff(u(t,x),t)-u(t,x)^m=0;
pde:=diff(u(t,x),t)-u^m=0;
pdsolve(pde); #Both work
pdsolve(pde2); #
## No difference with declare except for the compact display.
restart;
with(PDEtools,declare);
declare(u(t,x));
pde2:=diff(u(t,x),t)-u(t,x)^m=0;
pde:=diff(u(t,x),t)-u^m=0;
pdsolve(pde);
pdsolve(pde2);
################
This behavior is rather inconsistent within pdsolve, but certainly also with dsolve, where arguments cannot be left out in input.
I shall submit an SCR on the inconsistency.

Your ode (vgl) is linear and homogeneous, thus the question about a particular solution doesn't come up.

That Maple's dsolve cannot find the general solution, but returns a DESol structure is just a fact.
You could find two linearly independent solutions by numerical solution:
res1:=dsolve({vgl,y(0)=1,D(y)(0)=0},numeric);
res2:=dsolve({vgl,y(0)=0,D(y)(0)=1},numeric);

or you could find terms in a series expansion for those solutions:

dsolve({vgl,y(0)=1,D(y)(0)=0},y(x),series,order=10);
dsolve({vgl,y(0)=0,D(y)(0)=1},y(x),series,order=10);

After having executed everything till just before dsolve, try
indets({Lagran1, Lagran2, Lagran3, Lagran4, Lagran5},function);

You will see weird constructions like cos((q2(t))(t)), diff((q1(t))(t), t, t).

Secondly, in ini you should use D for the derivatives at a point. Thus q1'(0) should be written D(q1)(0). The second derivative at 0 as (D@@2)(q1)(0). ##For the second derivative, see the note at bottom.

#### The first time the "weird" construction comes up is here:
Eq53 := subs(q1=q1(t),q2=q2(t),q3=q3(t),q4=q4(t),q5=q5(t),v1=diff(q1(t),t),v2=diff(q2(t),t),v3=diff(q3(t),t),v4=diff(q4(t),t),v5=diff(q5(t),t), Eq15); #replace your colon by a semicolon to see that.
Did you actually mean to substitute into Eq51 instead?
As it happens Eq51 and Eq52 are both zero.
When I use Eq53 as defined by a substitution into Eq51 (which is 0) instead of into Eq15, then I get
indets({Lagran1, Lagran2, Lagran3, Lagran4, Lagran5},function);
where q5 doesn't appear at all, so causes the following error from dsolve:

Error, (in dsolve/numeric/process_input) indication of the unknown(s) of the problem {q1(t), q2(t), q3(t), q4(t), q5(t)} disagrees with input system {q1(t), q2(t), q3(t), q4(t)}

Notice that then Lagran5 is the equation 0=0. Since you don't have to give the functions as the second argument you could try:
dsolve({Lagran1, Lagran2, Lagran3, Lagran4, Lagran5, ini}, numeric, output=listprocedure);
then you get the unsurprising error

Error, (in dsolve/numeric/process_input) invalid specification of initial conditions, got {0 = 0, q5(0) = 0, (D(q5))(0) = 0}
### An additional problem:
Eq34 := subs(vq1=q1(t), .....)
You probably meant q1 not vq1.
## And one more: For a second order systen you cannot impose initial conditions on the second derivatives.

From a cursory look at the procedure sin (and from some experimentation) it seems to me that lexical ordering may be the explanation for the behavior, not some preference for minus over plus.

showstat(sin);

In line 40 of the procedure we find:

elif type(x,`+`) and `tools/sign`(x) = -1 then
  40     res := -sin(-x)

You may try
`tools/sign`(z-y);
and you get -1.

If you look at
showstat(`tools/sign`);

you will see that sign is being used in line 7 and
sign(z-y);
returns -1.
Now according to the help page for sign, "The sign function computes the sign of the leading coefficient of expr".
Furthermore it states:
"The leading coefficient of expr is determined with respect to the indeterminates given.  If none are given, the leading coefficient is taken with respect to all its indeterminates.  Note therefore that the leading coefficient is dependent on the order of the indeterminates which may vary from one Maple session to another, but not within a session." (my emphasis).

lcoeff(z-y); # I got -1

Finally in the help page for lcoeff you find that the order from indets is used:
indets(z-y);  # I got {y,z}
## Actually the page makes this more specific:
frontend(indets, [z-y], [{`*`, `+`, `::`, series, SDMPolynom, constant, undefined}]); # returns the same in this case

Begin by writing down the system in the syntax used by Maple.
Here is your first ode where for simplicity in input I have used the product rule of differentiation.

ode1:=diff(f(x),x,x)+1/5*diff(p1(x)*f(x),x)+2/5*diff(p2(x)*g(x)*f(x),x)=0;

Similarly write down ode2 and ode3. (I hope I got the typing OK: IT IS A BORE!, so next time you do it!).

ode2:=diff(h(x),x,x)*f(x)+2*diff(h(x),x)*diff(f(x),x)+1/5*p1(x)*f(x)*diff(h(x),x)+2/5*p2(x)*g(x)*f(x)*diff(h(x),x)-1/5*diff(p2(x)*f(x)^2,x)=0;
ode3:=diff(h(x),x)=1/5*p2(x)*f(x);

Use dsolve:

res:=dsolve({ode1,ode2,ode3},{f(x),g(x),h(x)});
  

The result is a sequence of two solutions. The first is pretty simple and rather trivial: f identically zero, h an arbitrary constant, and g can be any function.

In the second solution, f can be any function, and h and g are expressed in terms of f, p1, and p2.

Clearly a bug. (I submitted an SCR).

A workaround is to apply ln on both sides of the equation to get x*ln(x+exp(-1)) = 0 and then solve:

solve(x*ln(x+exp(-1)) = 0, x); #Two solutions found: 0, 1-exp(-1))
## Verification:
eval((x+exp(-1))^x,x=1-exp(-1));
eval((x+exp(-1))^x,x=0);
# Graphical illustration:
plot((x+exp(-1))^x -1, x=-exp(-1)..2);
## fsolve can find the two real roots easily:
fsolve((x+exp(-1))^x = 1, x);
fsolve((x+exp(-1))^x = 1, x=0.5);
## It seems that the original equation also has imaginary roots:
plots:-complexplot3d((x+exp(-1))^x - 1,x=-3-7*I..3+7*I,view=0..2);
fsolve((x+exp(-1))^x = 1, x=2+3*I,complex);
fsolve((x+exp(-1))^x = 1, x=2-3*I,complex);
fsolve((x+exp(-1))^x = 1, x=3+5*I,complex);
fsolve((x+exp(-1))^x = 1, x=3-5*I,complex);
####
Note. (Thanks to Axel Vogt).
Using a more careful reformulation of the original equation we can state that
the equation (x+exp(-1))^x = 1 is equivalent to x*ln(x+exp(-1)) = 2*Pi*k*I for some integer k.
Solving the latter we don't get an explicit formula, but a RootOf expression, which can be used to draw a graph of some of the infinitely many solutions.
k=0 is a special case, and is already handled above. Notice that we include 1-exp(-1) explicitly in the plot.
r:=solve(x*ln(x+exp(-1)) = 2*k*I*Pi, x);
plots:-complexplot([1-exp(-1),seq(r,k=-10..10)],style=point,symbol=solidcircle,symbolsize=20);

As Carl mentions, dsolve certainly handles bvp problems numerically. As an example I take your system.
I start by the corrected version of ODE mentioned in my answer about shooting.
After that I reconstruct the original system of higher order equations, not because it is necessary, but to illustrate that rewriting the original system as a first order system can be left to dsolve: the user doesn't have to do it.

ODE2 := {diff(f(eta), eta) = fp(eta), diff(fp(eta), eta) = fpp(eta), diff(fpp(eta), eta) = fppp(eta), diff(fppp(eta), eta) = -96.*fp(eta)*m(eta)^2+96.*mp(eta)*f(eta)*m(eta)-48.*mp(eta)*m(eta)*eta-2.*fppp(eta)*f(eta)+fppp(eta)*eta-48.*m(eta)^2+3.*fpp(eta), diff(g(eta), eta) = gp(eta), diff(gp(eta), eta) = 2.*g(eta)*fp(eta)-2.*gp(eta)*f(eta)+gp(eta)*eta-2.*n(eta)*mp(eta)+2.*m(eta)*np(eta)+2.*g(eta), diff(m(eta), eta) = mp(eta), diff(mp(eta), eta) = 12.*m(eta)*fp(eta)-12.*mp(eta)*f(eta)+6.*mp(eta)*eta+6.*m(eta), diff(n(eta), eta) = np(eta), diff(np(eta), eta) = 48.*m(eta)*gp(eta)-12.*np(eta)*f(eta)+6.*np(eta)*eta+12.*n(eta)};
##
## Now take the extra equations you introduced as a user in order to write the system as a first order system:
##
peqs:={diff(f(eta), eta) = fp(eta), diff(fp(eta), eta) = fpp(eta), diff(fpp(eta), eta) = fppp(eta),diff(g(eta), eta) = gp(eta),diff(m(eta), eta) = mp(eta),diff(n(eta), eta) = np(eta)};
##
sys1:=ODE2 minus peqs;
sbs:=evalindets(peqs,`=`,x->rhs(x)=lhs(x));
## Your original system was the following (or an equivalent):
sys:=eval[recurse](sys1,sbs);
## Now reconstruct the boundary conditions from IC and BC:
##
BCS:= {f(0) = 0, g(0) = 1, m(0) = 0, n(0) = 0, D(f)(0) = 0, f(blt) = .5, D(f)(blt) = 0, D(g)(blt) = 0, m(blt) = 1, n(blt) = 1};
##
## Now solve the bvp problem:
res:=dsolve(sys union BCS,numeric);
## Example plots:
plots:-odeplot(res,[seq([eta,diff(f(eta),[eta$i])],i=0..3)],0..blt);
## Notice that you won't get away with plotting over an interval larger than 0..blt due to the method used. Normally that is not what people would want to do anyway.

Replace 'recursive' by 'factor'.  (recursive is default).

In Maple 2016 I don't get the error you report, but I get another one:
Error, (in isolate) unable to isolate YP[11]

This error comes from `dsolve/numeric/bvp/convertsys` when it does this:
`dsolve/numeric/bvp/convertsys`(eval({Eq1, Eq2, Eq3, Eq4, Eq5, Eq6},Pr=L[1]),{},[F,f,phi,phi1,theta,theta1],eta,l);
# There is no problem with DEtools[convertsys] which is not used by `dsolve/numeric/bvp`:
DEtools[convertsys](eval({Eq1, Eq2, Eq3, Eq4, Eq5, Eq6},Pr=L[1]),{},[F,f,phi,phi1,theta,theta1],eta);

As Tom Leslie found out, this is a case of garbage in, garbage out.