Clearly a bug. And an old one at that. It appears in Maple 8 and Maple 12 as well as in Maple 2017.2.

The problem may be located in line 32 of `singular/singular` which reads:
assume(new_Z,'posint');

## I shall submit an SCR (aka a bug report).

Take a look at the help page

?Initially Known Names

You will see that gamma(0), gamma(1), ..., gamma(5), ... are known as Stieltjes constants. Furthermore gamma = gamma(0) is Euler's constant.
That is the problem, I'm sure. The unprotecting you are doing doesn't change that unless you follow it up by assigning gamma to some other name e..g. gamma:=g;
But of course that is silly: just use g to begin with instead of gamma (or h as Tom Leslie is using).

An alternative in recent versions of Maple is to place right after restart the local declaration:
local gamma;
## Omit the unprotect gamma part: it serves no purpose.
#### To see the numerical values of gamma(n) for n=0,1,2... you can try
seq(evalf(gamma(n)),n=0..5);
## To see the code for this, try
showstat(`evalf/gamma`);

The mathematical constant in question is in Maple spelled Pi, not pi. The latter is just another symbol with no meaning, but it prints the same as Pi does.
cos(Pi/2)^2;

will obviously just evaluate to 0.

You can use a try statement as in
 

Q1 := []; 
for jj to 100 do 
  simplfloat := rand(-1.0 .. 1.0); 
  h := simplfloat(); 
  a := abs(amin+(amax-amin)*h); 
  b := abs(bmin+(bmax-bmin)*h); 
  try 
    solut := dsolve([eq, cis], numeric, maxfun = 0, output = procedurelist); 
    sd := i-> abs(rhs(solut(i)[2])*rhs(solut(i)[2])); 
    eng := Simpson(sd, 5, 6, 10) 
  catch "cannot evaluate the solution": 
    printf(StringTools:-FormatMessage(lastexception[2 .. -1]));
    printf("\n"); 
    eng := FAIL 
  end try; 
  if eng = FAIL then next else Q1 := [op(Q1), [jj, a, b, h, eng]] end if 
end do:

But in my attempts with the loop going to 100 I never get any successes. So the question is: Must a singularity occur before the end of the interval?
Since what seems to happen is that x(t) decreases to zero thus creating a problem with 1/x(t), we may rephrase the question to: Must x(t) decrease to zero under the given conditions?
### Actually, I think the answer is yes, it must.
Rewriting the ode as a first order system using the usual diff(x(t),t) = y(t), we get:
 

sys:=[diff(x(t), t) = y(t), diff(y(t), t) = -(b*y(t)+sin(x(t))*y(t)+2/x(t))/a];

A small amount of phase plane analysis making use of the vertical isocline y = 0 and the horizontal isocline y = -2/x/(b+sin(x)) and the fact that a and b are positive with b > 1 will show that if the solution starts at (x,y) = (0.25, 0) (i.e. on the vertical isocline) then it will not be able to cross the horizontal isocline, which is situated in the 4th quadrant. Thus the only possibility is that x(t) -> 0 and y(t) -> -infinity as t increases to the right end of the (finite) interval of definition of the solution.
Illustration here for a = .1828083187 and b = 22.45274868.
The red line is the vertical isocline, the blue the horizontal isocline. The trajectory must stay between these thoughout its life (for t >=0 ).

Using DEplot:

Produced by:
 

plots:-display(
   DEtools:-DEplot(eval(sys,{a = .1828083187, b = 22.45274868}),[x(t),y(t)],t=0..0.38,[[x(0)=0.25,y(0)=0]],linecolor=maroon,color=gray),
   plot(-2/x/(22.45274868+sin(x)),x=0..0.25,color=blue),
   plot(0,x=0..0.25,color=red,thickness=3),
   view=-2.5..0
);

Let us consider the simplified example sol:=1/x/sin(x) to understand what is going on.
If you look at the code for algsubs:
showstat(algsubs);
you will find that you are going to `algsubs/algsubs`  with the call

`algsubs/algsubs`(1/x/sin(x),x,t,[x],remainder);

Now look at the code for `algsubs/algsubs`:
showstat(`algsubs/algsubs`); #Very short
You will see that you are going to `algsubs/recurse` with the call

`algsubs/recurse`(1/x/sin(x),x,t,[x],remainder);

Now looking at the code for `algsubs/recurse`:
showstat(`algsubs/recurse`);
you will that since type(1/x/sin(x),`*`) = true, you will be in line 4 with
map(`algsubs/recurse`, 1/x/sin(x), x, t,[x], remainder);
Thus you will handle both of these
`algsubs/recurse`, 1/x, x, t,[x], remainder);
`algsubs/recurse`, 1/sin(x), x, t,[x], remainder);
The first one will be handled by `algsubs/recurse` in line 3 since
type(1/x,'ratpoly'(anything,[x])) = true
so 1/x is returned for that factor.
As for 1/sin(x) that will be handled in line 7 of `algsubs/recurse` with the call

map(`algsubs/algsubs`,1/sin(x),x,t,[x],remainder);
resulting in 1/sin(t).
Thus altogether 1/x/sin(t).
So the crucial line is line 3 in `algsubs/recurse`:
elif type(f,('ratpoly')(anything,vars)) then f

As Kitonum says, this is not a bug. In my view because it is so blatant: It must have been intended.
As Kitonum also points out, there is this note in the help page:

"Note: The algsubs command currently works only with integer exponents."

I find that rather obscure: Integer exponents of what? And isn't -1 an integer?

## My own note: I never use algsubs, only subs and eval.

In the following I assume that the syntax errors (missing multiplication signs) have been corrected.
eq2 and eq3 are linear homogeneous odes and are subject to homogeneous boundary conditions.
So clearly psi = w = 0 (identically) is a solution. Is it the only one?
## I have looked at that question now. It is the only solution.
I have placed the proof at the end of this.
##
But let us look at that one.
So we have:
 

restart;
eq1 := (diff(psi(x), x))^2+(diff(u(x), x)+(8*(1/2))*(diff(w(x), x))^2)*(diff(psi(x), x))^2+3*(diff(w(x), x, x))+5*(diff(w(x), x, x))*(diff(psi(x), x))-7*(diff(u(x), x, x, x)+(8*(1/2))*(diff(w(x), x, x))^2+(3/2)*(diff(w(x), x, x, x))*(diff(w(x), x)))+3 = p;
eq2 := (51-31)*(diff(psi(x), x, x))+(52-2)*(diff(w(x), x, x, x))+8*(diff(psi(x), x, x, x, x))-7*(diff(w(x), x)-psi(x)) = 0;
eq3 := 4*(diff(w(x), x, x)-(diff(psi(x), x)))+(23+11)*(diff(psi(x), x, x, x))+(14+12)*(diff(w(x), x, x, x, x)) = 0;
bcs:={psi(0) = 0, psi(1) = 0, u(0) = 0, u(1) = 0, w(0) = 0, w(1) = 0, D(u)(1) = 0, D(psi)(0) = 0, D(psi)(1) = 0, D(w)(0) = 0, D(w)(1) = 0};
bcs1,bcs23:=selectremove(has,bcs,u);
eq10:=eval(eq1,{psi=0,w=0});
sol0:=dsolve({eq10} union bcs1);
plots:-animate(plot,[rhs(sol0),x=0..1],p=-50..50);

So we see that for each p there is a solution for u, psi, w where psi = w = 0, and u<>0 unless p=3:
 

solve(identity(rhs(sol0),x),p);
eval(sol0,p=3);

So you may ask, what happens numerically? Well, begin by plotting w and psi in my previous answer, where I only plotted u.
You will find that psi = w = 0.

### I place here the proof that eq2, eq3 with bcs23 only has the trivial solution psi = 0, w = 0.
 

restart;
eq1 := (diff(psi(x), x))^2+(diff(u(x), x)+(8*(1/2))*(diff(w(x), x))^2)*(diff(psi(x), x))^2+3*(diff(w(x), x, x))+5*(diff(w(x), x, x))*(diff(psi(x), x))-7*(diff(u(x), x, x, x)+(8*(1/2))*(diff(w(x), x, x))^2+(3/2)*(diff(w(x), x, x, x))*(diff(w(x), x)))+3 = p;
eq2 := (51-31)*(diff(psi(x), x, x))+(52-2)*(diff(w(x), x, x, x))+8*(diff(psi(x), x, x, x, x))-7*(diff(w(x), x)-psi(x)) = 0;
eq3 := 4*(diff(w(x), x, x)-(diff(psi(x), x)))+(23+11)*(diff(psi(x), x, x, x))+(14+12)*(diff(w(x), x, x, x, x)) = 0;
bcs:={psi(0) = 0, psi(1) = 0, u(0) = 0, u(1) = 0, w(0) = 0, w(1) = 0, D(u)(1) = 0, D(psi)(0) = 0, D(psi)(1) = 0, D(w)(0) = 0, D(w)(1) = 0};
bcs1,bcs23:=selectremove(has,bcs,u);
### Now considering only eq2 and eq3 with bcs23
sys23:={eq2,eq3};
bcs23;
## Splitting in conditions at x = 0 and x = 1:
bcs23_0,bcs23_1:=selectremove(hastype,bcs23,anyfunc(0));
## Defining initial conditions at x = 0 using names for unknown values:
ics23:=bcs23_0 union {(D@@2)(w)(0)=w2,(D@@3)(w)(0)=w3, (D@@2)(psi)(0)=p2,(D@@3)(psi)(0)=p3};
nops(ics23); # 8 as there should be
## We now have a standard initial value problem at 0 for sys23 with ics23.
## This has a unique solution for any w2, w3, p2, p3:
sol:=dsolve(sys23 union ics23): # This fills a lot, therefore the colon.
evalf(sol);# Just to see a smaller output, but will NOT be used
## Now imposing the conditions at x = 1 on the solution:
EQ1:=eval[recurse](sol,{x=1} union bcs23_1): 
evalf(EQ1); # Just having a look at a smaller version.Will NOT be used. 
## Some of the conditions at x=1 involve the first derivatives:
EQ2:=eval[recurse](convert(diff(sol,x),D),{x=1} union bcs23_1):
evalf(EQ2);# Just having a look at a smaller version.Will NOT be used.
indets({EQ1,EQ2},name);
## Now rewriting the linear system EQ1, EQ2 as a matrix equation: 
M,zz:=LinearAlgebra:-GenerateMatrix(EQ1 union EQ2,[p2,p3,w2,w3]):
zz;
evalf(M); # Just having a look at a smaller version.Will NOT be used.
## The determinant of M:
DT:=LinearAlgebra:-Determinant(M):
evalf(DT); # Apparently different from zero
evalf[50](DT); # Confirmed
is(DT>0); # Confirmed
## Thus the system M.< p2,p3,w2,w3 > = <0,0,0,0> has a unique solution.
## Since < p2,p3,w2,w3 > = <0,0,0,0> is a solution, it is the one.

My guess is that it is a simple bug. Try this:

with(ScientificConstants);
debug(Element);
GetValue(Element('H',density(gas)));
GetValue(Element('Si',density));
showstat(Element);

You will see that in line 47 of the procedure Element:
showstat(Element, 47);
ScientificConstants:-numTabElem[an]:-props[yp]:-value(`if`(withPs,params,NULL))
in the case of Hydrogen produces the number 0.088 because withPs = true, whereas for Silicium NULL is produced because withPs=false.
So my guess is that for gases the numbers produced at this stage are actually meant to be in kg/m^3, but are incorrectly treated as if in g/cm^3 and thus converted to kg/m^3 later (by GetValue in its line 42).
A closer look at line 42 in the procedure GetValue shows that in
convert(valueX,'system',ScientificConstants:-numTabElem[op(x)[1]]:-props[op(0,op(x)[2])]:-units,sys)
this part:
ScientificConstants:-numTabElem[op(x)[1]]:-props[op(0,op(x)[2])]:-units
evaluates to kg/L and sys evaluates to SI. valueX has the value 0.088.
Now kg/L is these days exactly the same as g/cm^3.
The numbers for gases must have been entered in SI units and are here are mistakenly treated as being in kg/L.

I shall submit an SCR with a link to this page.
PS. In preparing for the SCR I discovered that the density in SI units for Hydrogen in Maple 12 and Maple 15 is given as 0.000088. This changed to 88.0 in Maple 16 and all later versions.

The error message tells the story: You have a system of total order 4+4+3 = 11 and you have an unassigned parameter p.
If your intention is to determine p at the same time you solve, you must have one more boundary condition. That makes it a total of 12.

You are using the indexed variables a[i] and b[i], but also the unrelated (your intention, I assume) a and b.
That is not a good idea. As soon as you do e.g.
a[1] := 0;
a table by the name 'a' is implicitly created.
Thus it is potentially catastrophic to define as you do
a[-1] := -12*mu/(a+b);
Also a and b appear in your pde.
The solution is simple: Either change the names a and b to e.g. aa and bb (or whatever) in the pde and in
a[-1] := -12*mu/(aa+bb);
or don't use indexed variables. Use other names for a[0], a[-1], etc.
### To see the catastrophy in a simple session try this:
 

restart;
a[-1] := -12*mu/(a+b);
eval(a); #table with name a and index -1 refers to the whole table itself
simplify(a); # Not surprising!

Yes, this is a bug. You are doing nothing wrong.
The two cases are treated by separation of variables: Try
 

pdsolve(PDE1);
pdsolve(PDE2);

Both of these results are correct, just different: _c[1] in the first corresponds to c[1] - a in the second. _c[1] is an arbitrary constant at this stage.
You could go further and ask for solutions of the odes:

pdsolve(PDE1,build);
pdsolve(PDE2,build);

The solutions are on different forms for sure, but they are still solutions of their respective pdes.
So whatever goes wrong happens later.
You can if you wish try your commands after the command
debug(pdsolve);
##
I shall submit a bug report with a reference to this question.

Your problem is that this command (solving for the highest derivatives):
 

solve(newsys2,{diff(g1(y), y$5),diff(g2(y), y$4),diff(g3(y), y$6),diff(g4(y), y, y),diff(g5(y), y, y)});

has no solution.
One of the equations (in the numbering on my machine it is newsys2[4]) contains none of the highest derivatives, i.e. none of
diff(g1(y), y$5),diff(g2(y), y$4),diff(g3(y), y$6),diff(g4(y), y, y),diff(g5(y), y, y).
There is a simple solution: Differentiate that equation w.r.t. y and use that instead of newsys2[4].
Thus your new system with the highest derivatives solved for is:

SYS:=solve({newsys2[1],newsys2[2],newsys2[3], diff(newsys2[4], y), newsys2[5]}, {diff(g1(y), y$5),diff(g2(y), y$4),diff(g3(y), y$6),diff(g4(y), y, y),diff(g5(y), y, y)});

Remember though that since diff(newsys2[4],y)=0 we have that  newsys2[4] = constant, thus demanding that newsys2[4] be satisfied at one of the boundary points must be added. Thus either add this boundary condition:
 

eval[recurse](convert(newsys2[4],D),{y=0} union bcs3);

or this one
 

eval[recurse](convert(newsys2[4],D),{y=1} union bcs3);

## A comment on the boundary conditions.
You have a system that when converted to a first order system has 19 equations. You have an undetermined constant omega3, so you need 19+1 = 20 boundary conditions.
bcs3 has 18 conditions and one of the two mentioned above (from newsys2[4] at 0 or 1) makes a total of 19. You need one more condition. The extra condition can involve the following quantities:
 

N:=[5,4,6,2,2];
extra_bcs := {seq(seq(seq((D@@i)(cat(g, j))(a), i = 0 .. N[j]-1),j = 1 .. nops(N)), a = 0 .. 1)};

Since neither the conditions in bcs3 nor in the newly added one are simple i.e. of the form (D@@i)(gj)(a) = constant (a=0 or 1) it doesn't help much to define extra_bcs with an added minus bcs3.

It appears that you are using 2D math input and have chosen to define f as a function.

You might just as well have defined your expression as just that: an expression in x, like this:
f1:=2*sin(0.05*Pi*x-Pi*0.5)+2;

If you don't give a range, plot will choose one for you, sometimes -10..10 under some circumstances -2*Pi..2*Pi.
Try e.g.
plot(cos);
In your case do
plot(f, 0..40);
Alternatively:
plot(f(x), x = 0..40);
###
In the expression case (using f1) just do
plot(f1, x = 0..40);
###
The rewriting you are talking about must refer to the default treatment of Pi in products with a float, but that is irrelevant for your present plotting problem.
Try
Pi*0.5;
If you get 1.570796327 you have the default setting of floatPi=true in kernelopts.
That can be set to false, as I myself have done.
kernelopts(floatPi=false);
You can put that statement into your maple.ini file, if you wish.

You say that exp(alpha1*(Dp*alpha1*t-lh+x))  is infinite in one interval. That could only be the case if alpha1*(Dp*alpha1*t-lh+x) is infinite in one interval. But that quantity is finite for all t >=0 (and indeed all t, but here only positive t are interesting).

(You have a thetac and a Dc in your assumptions, but they are not in the expression.)

Basically what you have is:
 

eq1:=exp(-a*sqrt(s))/(sqrt(s)+b);
inttrans[invlaplace](eq1,s,t) assuming a>0,b>0;

which returns
exp(-(1/4)*a^2/t)/sqrt(Pi*t)-b*exp(b*(b*t+a))*erfc((1/2)*(2*b*t+a)/sqrt(t))

What you have in your definition of A and B are sequences. A:=(-4,5,-2) defines a sequence. You don't even need the parenthesis.
What you need to do is to define vectors e.g. by using <... > as in:
 

A,B:=<-4,5,-2>,<8,2,7>;
A-B;

There are several syntactical errors:
1. While odes is a sequence, ics is a set, so you must change the dsolve command, e.g. to dsolve({odes} union ics, .... ).
2. Did you intend to set Digits:=50? in that case change digits:=50 to that.
3. What was the intended meaning of the 4th ode:
 

ode4:=diff(r^4*exp(-v(r))*sqrt(1-2*G*M(r)/(r*c^2))*(diff(1-omega(r)/Omega, r)), r)+4*r^3*(1-omega(r)/Omega, r)*(diff(exp(-v(r))*sqrt(1-2*G*M(r)/(r*c^2)), r)) = 0;
## This part in particular:
op(5,lhs(ode4));
nops(%); # A product of 4 factors

The product of the first 3 factors is 4*r^3*(1-omega(r)/Omega, r). The last factor doesn't make any sense in this context.