In Maple 17 you can recapture the old result by doing this:
restart;
T := (4*(a-(1/2)*c-(1/2)*b)^2*cos((2/3)*xi*sqrt(a^2+(-c-b)*a+c^2+b^2-b*c)/r)+3*(-c+b)^2)/(4*b^2+(-4*a-4*c)*b+4*a^2+4*c^2-4*a*c);
f:=unapply(T,a,b,c);
g:=t->f(a*t,b*t,c*t);
series(g(t),t=0,5);
eval(convert(%,polynom),t=1);
res:=collect(simplify(%),[a,b,c],distributed);
####################################
This behavior seems to have begun in Maple 16 and is not a bug in mtaylor, but rather a weakness (or bug) in taylor itself. taylor is builtin, whereas mtaylor is not.
####################################
A very simple version:
restart;
V:=cos(sqrt(x^2));
taylor(V,x,5); #Works in Maple 15 but not in Maple 16 or 17
series(V,x,5);
simplify(%);
#############
I shall submit an SCR of the latter simple example.

restart;
with(inttrans):
T:=tanh(s);
T1:=normal(1-convert(T,exp));
#Can be written:
T2:=2*exp(-2*s)/(1+exp(-2*s));
#Can be written as an infinite series:
T3:=2*Sum((-1)^n*exp(-2*s*(n+1)),n=0..infinity);
value(T3); #Check
#Applying termwise invlaplace
f3:=evalindets(T3,specfunc(anything,exp),u->simplify(invlaplace(u,s,t),assume=[n>=0]));
#What is going on there was
invlaplace(exp(-2*s*(n+1)), s, t)  assuming n>=0;
#So the result should be
f:=invlaplace(1,s,t)-f3;
#Check
laplace(f,t,s);
value(%);
simplify(%-convert(T,exp));

I have changed both of your procedures some. They both return a vector:

restart;
procVec1 := proc (U) local m,N,i;
m := nops(U);
N := Vector(m-1);
for i to m-1 do
    N[i] := unapply(piecewise(U[i] <= u and u < U[i+1], 1),u);
end do;
N
end proc;
procVec2 := proc (U) local m;
m := nops(U);
Vector(m-1,i->unapply(piecewise(U[i] <= u and u < U[i+1], 1),u))
end proc;
U:=[0,1,2,3,4];
procVec1(U);
procVec2(U);

In the real domain (which we are dealing with) the range of arccos is [0,Pi].
Thus there are no solutions to
{x > 0, y > 0, 30 = sqrt(x^2+y^2), a = arccos(abs(x)/sqrt(x^2+y^2))}
when a < 0.
So clearly the answer from Maple 17 is wrong.
The range of the restriction of arccos to [0,1] is [0,Pi/2], thus there are in fact only solutions for a in [0,Pi/2].
The solution is x = 30*cos(a), y = 30*sin(a).
#Check:
eval({x > 0, y > 0, 30 = sqrt(x^2+y^2), a = arccos(abs(x)/sqrt(x^2+y^2))},{x=30*cos(a),y=30*sin(a)});
simplify(%) assuming a>=0,a<=Pi/2;

I tried shooting and, as it seems, successfully.
The idea is to replace the boundary conditions with the initial conditions
ics:=f(h1) = (1/2)*F, D(f)(h1) = f1, (D@@2)(f)(h1)=f2,(D@@3)(f)(h1)=f3,
      theta(h1)=-g*th1, D(theta)(h1)=th1 ;
 and then use dsolve with parameters f1, f2, f3, th1.
Then define four procedures that must return 0 for the boundary conditions to be satisfied.

MaplePrimes13-08-17s.mw

I may add that after having solved by the shooting method with the initialized res you can give res to the boundary problem as an approximate solution:

Bsys:=[Eq1,Eq2,bcs];
resB:=dsolve(eval(Bsys,param),numeric,approxsoln=res);
plots:-odeplot(resB,[[y,f(y)],[y,theta(y)]],H2..H1);


It is comforting but hardly surprising that the result comes right away.

The 3rd and fourth of your boundary conditions are cubic polynomials in f''. Thus in solving for f'' you may have 1 or 3 real solutions depending on f'.
Take the third boundary condition:
bcs3:=D(f)(h1)+beta*(alpha*(D@@2)(f)(h1)-beta2*(D@@2)(f)(h1)^3) = -1;
#There is exactly one real root for f''(h1) iff the discriminant of the following polynomial in f2=(D@D)(f)(h1) is nonnegative:
pol:=eval((lhs-rhs)(bcs3),{(D@@2)(f)(h1)=f2,D(f)(h1)=f1});
dsc:=discrim(pol,f2);
eval(dsc,{beta = 0.5, beta2 = 0.1,alpha=2});
#So you see that the number of real roots depend on f1 = D(f)(h1):
sol:=eval([solve](pol,f2),{beta = 0.5, beta2 = 0.1,alpha=2});
evalf(eval(sol,f1=1));
simplify(fnormal(evalf(eval(sol,f1=0))));
#The same applies to the fourth boundary condition.
It is probably worth trying to solve for f2 and pick 1 of the 3 roots before asking dsolve to do anything.
########
In the help page dsolve/numeric/BVP it says
"The solution method is capable of handling both linear and nonlinear BVPs with fixed, periodic, and even nonlinear boundary conditions. For nonlinear boundary conditions, however, it may be necessary to provide an initial solution profile that approximately satisfies the boundary conditions"
I got the error message:
"Error, (in dsolve/numeric/bvp) cannot determine a suitable initial profile, please specify an approximate initial solution"
when I tried picking the first solution f2 at both ends. Even when omitting theta entirely.
##############################################################
First time around I only looked at the output from Eq1. However, writing Diff instead of diff like this:
Eq1:=Diff(alpha*diff(f(y),y$2)-M^2*f(y)-beta2*(diff(f(y),y$2))^3,y$2)=0;
revealed to me (finally) that actually your first equation could be written
Eq1a:=alpha*(diff(f(y), y, y))-M^2*f(y)-beta2*(diff(f(y), y, y))^3=c1*y+c2;
where c1 and c2 are undetermined constants. Thus you have an ode cubic in f''. Could you proceed from there?

Try using 1D-math (aka Maple Input).
My default is worksheet interface and Maple Input (set via the menu Tools/Options etc.), but I have the same problem with my keyboard if I ask for 2D-math input either by using the menu 'Insert' or by doing CTRL R.
I have no problem with 1D input.

And like you I don't have the problem in Maple 16.

A workaround which works for me is to use the ascii code 125 for }. I hold down the Alt-key and enter 125 on the numeric keyboard to my right. That gives me }. (The other { is 123).

I shall report the problem by submitting an SCR.

Here is a somewhat different approach. The first uses dsolve directly, the other shows some intermediate steps:

restart;
eq10:=diff(h(t), t)*h(t) = 1/8*r^2*(hmax-h(t))*rho*g*sin(alpha)/eta;
res1:=rhs(dsolve({eq10,h(0)=0}));
#Second approach:
dsolve({eq10,h(0)=0},useInt,implicit);
IntegrationTools:-Expand(%);
value(%) assuming hmax>0,h(t)<hmax,h(t)>0;
res2:=solve(%,h(t));
#The results are indeed the same:
evalindets(res1-res2,specfunc(anything,exp),expand);

In your loop you are defining functions (procedures). At the time the function is defined nothing is evaluated. Thus i is not evaluated.
To force evaluation use unapply:

restart;
for i from 0 to 7 do
  gg[i]:=unapply(t^i,t)
end do;
gg[3](u);
gg[0](u);

Converting to rationals and using diff helps:
pde:=(D[1](u))(t, x)-(1/2*(x-x^2))*(D[2, 2](u))(t, x)-(-5.02*x+4.51+.1*x*(1-x))*(D[2](u))(t, x) = 0;
pde2:=convert(convert(pde,rational),diff);
pdsolve(pde2);
pdsolve(pde2,build);
#Numerical solution
res:=pdsolve(pde2,{u(t,0)=.1,u(t,1)=0,u(0,x)=0},numeric);
#Another example. Notice in particular that boundary conditions are given at x = 0.1 and x = 1:
res:=pdsolve(pde2,{u(t,0.1)=0,u(t,1)=0,u(0,x)=x^2},numeric,spacestep=.01);
res:-animate(t=2);

This is undoubtedly a bug in simplify.
restart;
assume(k>0,hbar>0,mu>0);
omega := hbar*sqrt(k/mu);
ycon := sqrt(sqrt(k*mu/hbar^2));
y := x-> ycon*x;
eta := (v, x)->sqrt(ycon/(2^v*factorial(v)*sqrt(Pi)))*HermiteH(v, y(x))*exp(-(1/2)*y(x)^2);
U:=(-hbar^2*(diff(diff(eta(0, x), x), x))/(2*mu)+(1/2)*k*x^2*eta(0, x))/(omega*eta(0, x));
#In order to get closer to the problem we may make it somewhat simpler while keeping the bug:
V:=eval(U,{HermiteH=1,hbar=1,x=1,exp=1});
simplify(numer(V))/simplify(denom(V)); #OK
normal(V); #OK
simplify(%); #OK
factor(V); #OK
simplify(%); #OK
simplify(V); #NOT OK
simplify(eval(V,{k=2})); #OK

I already posted an answer to your first version:

http://www.mapleprimes.com/questions/150381-Closed-Piecewise-Function-In-3dimensional

You could use allvalues:

R:=indets(newpar,RootOf);
nops(R);
#The first of the roots in R[2] as an example:
allvalues(R[2])[1];

Maybe like this:

D124:=max(-x + y + z, x - y + z, x + y - z) = 1;
D134:=max(-x - y - z, -x + y + z, x - y + z) = 1;
D123:=max(-x - y - z, -x + y + z, x + y - z) = 1;
D234:=max(-x - y - z, x - y + z, x + y - z) = 1;
t:=piecewise(x>0 and y>0 and z>0,D124,x<0 and y<0 and z>0,D134,x<0 and y>0 and z<0,D123,x>0 and y<0 and z<0,D234);
plots:-implicitplot3d(t,x=-1..1,y=-1..1,z=-1..1);

You could do something like this:

N:=6:
if `and`(seq(thing[i],i=1..N)) then ... end if;