In the case with 8 unknowns Maple couldn't find any solution and issues a warning.

In the case with 6 unknowns you have 8 equations and only 6 unknowns. solve responds with NULL, which means no solutions.

So no reason for being surprised (if you are).

It has nothing to do with linalg although that package is very much outdated as Kitonum points out.

### Correction: It is not sys that is solved in the first solve command but a set of 6 equations (not being subject to the substitution used in defining sys), which is indeed linear in the chosen 6 variables. The coefficient matrix A, however, has rank only 5 while the augmented matrix <A|b> has rank 6, thus that linear system has no solution:
 

sys6:=[S__Fx = 0, S__Fz = 0, S__Myl = 0, S__Mys = 0, S__Mysi = 0, S__Mzle = 0];
nops(sys6); # 6
vars:=[R__xl, R__xs, R__xsi, R__zl, R__zs, R__zsi];
A,b:=LinearAlgebra:-GenerateMatrix(sys6,vars);
LinearAlgebra:-Rank(A); #5
LinearAlgebra:-Rank(<A|b>); #6

This is not a bug and there is no difference between Maple 2019.2 and Maple 2020.0 in the handling if the settings of floatPi are the same in the two versions.
What does the query kernelopts(floatPi); respond in your version of Maple 2019 and also in Maple 2020?

In my maple.ini file I have:
kernelopts(floatPi = false);

Notice that the default in recent versions after the (silly) idea was introduced is kernelopts(floatPi = true).

I never understood why the constant Pi should be treated any different from other constants like gamma, exp(1) or sqrt(2).
With kernelopts(floatPi = true). you get that when Digits = 10,

Pi/2-1e-12  =  1.570796327 in Maple 2019.2 as well as in Maple 2020.
Therefore we get in both versions
cos(Pi/2-1e-12)  =  -2.051033808*10^(-10).

If on the other hand kernelopts(floatPi = false) we get
Pi/2-1e-12  =  1/2*Pi - 1.*10^(-12)  in both versions, 2020 and 2019.2
and
cos (Pi/2 -1e-12) = 1.000000000*10^(-12) in both versions.

If you get anything different with Digits = 10 I should like to know.

#####

In Maple 18 the notion of of floatPi doesn't exist.
In Maple 2015 it was introduced (silently, not mentioned in the help for kernelopts). The default was kernelopts(floatPi=false).

The real introduction came with Maple 2016 and the default was kernelopts( floatPi=true).
It has remained that in 2017-2020.

####
To avoid an error in versions older than Maple 2015 and on the same computer what I have in my maple.ini file (besides other stuff) is
try kernelopts(floatPi=false) catch: end try;

Taking your desired answer and expanding it we have your starting point.

Use collect with extra argument factor:
 

E:=(2*h^2*sum(x[i]^2,i=1..4))*alpha+(2*h^2*sum(x[i]^2*y[i],i=1..4))*beta+2*h*sum(x[i]^2,i=1..4)-2*h*sum(x[i]*x[i+1],i=1..4);
expand(E);
collect(%,[alpha,beta],factor);

Here is a made up example of a 2x2 matrix ode:

restart;
U:=Matrix(2,2,symbol=u);
Ut:=apply~(U,t);
QL:=Matrix([[sin(t),cos(t)],[t,1]]);
RL:=Matrix([[t^2,ln(1+t)],[exp(-t),exp(-2*t)]]);
eq:=diff(Ut,t)=-QL.Ut-Ut.QL^%H-RL.Ut.RL^%H ;
sys:=Equate((lhs,rhs)(eq)) assuming real,t>-1;
init:=eval(Ut,t=0)=Matrix(2,2,symbol=u0);
inits:=Equate((lhs,rhs)(init));
res:=dsolve({sys[],inits[]},numeric,parameters=rhs~(inits));
res(parameters=[1,-2,7,0.1]);
res(5);
subs(res(5),Ut);
plots:-odeplot(res,[t,u[1,2](t)],t=0..5);

############
I chose to solve numerically, the more realstic approach in general.

You can change the initial conditions by just setting the parameters as in res(parameters=[1,-2,7,0.1]); above.

If initial conditions are fixed then omit the parameters option.

Replace the if statements with a nested piecewise.

piecewise returns unevaluated if x,y,z are not concrete.
 

azimuth:=piecewise(x>0,arctan(y/x),x<0,arctan(y/x)+Pi,x=0,piecewise(y>0,Pi/2,y=0,0,y<0,-Pi/2));

In Maple 2020 convert/Vector has been changed some:
Try
showstat(`convert/Vector`);

in earlier versions and compare with Maple 2020.

Immediately noticable is that the number of lines has basically doubled.

Let xxx be anything else than NULL, column or row.

In previous versions convert([1,2,3], xxx) ends up returning Vector[xxx]([1,2,3]), which results in the same error as the following does:

Vector[xxx]([1,2,3]);

In Maple 2020 the code for `convert/Vector`has in line 10 the following statement:
R := `if`(dir = NULL,Vector['column'](V),~Vector[dir](V));

which causes the problem.
Try in a new session:
 

restart;
dir:=xxx:
V:=[1,2,3];
R := `if`(dir = NULL,Vector['column'](V),~Vector[dir](V));

You will loose the kernel.

## Added later:  If you remove the "coercer" ~ from ~Vector[dir](V) above, you don't loose the kernel, but just get the appropriate error:

Error, (in Vector[xxx]) row or column expected for Vector index

This simple use of ~Vector causes a loss of kernel connection (in any Maple version from Maple 16 and up:
~Vector[xxx]([1,2,3]);

causes loss of kernel connection.

Surely you won't get anything useful out of trying to find an exact solution.

Using the values of R at eta = 0 you could get a series approximation that would (should) resemble the solution in a right-sided neighborhood of 0:
 

### For convenience in finding ics I use operator output in dsolve:
Digits:=15:
R := dsolve({bc, eq1, eq2}, numeric, output = operator,abserr=1e-12);
ics := R(0)[2 .. -1][];
S:=dsolve({eq1,eq2,ics},fcns,type=series,order=10);
Sp:=convert(S,polynom);
###
odeplot(R, [eta, f(eta) - eval(f(eta), Sp)], 0 .. 1); # OK on that short interval

######################################

Added later in the day:
An even better idea could be to solve the initial value problem using the ics we found from solving the bvp. Then find the differences. It should be noticed that the algorithms for solving bvp and ivp are quite different.

R_init:=dsolve({ics, eq1, eq2}, numeric, abserr=1e-12,relerr=1e-10,output=operator);
F,F1,F2,G,G1,G2:=op(rhs~(R_init[2..-1]));
odeplot(R,[eta,f(eta)-F(eta)],0..N);
odeplot(R,[eta,g(eta)-G(eta)],0..N);
odeplot(R,[eta,diff(f(eta),eta)-F1(eta)],0..N);
### Similarly with the other derivatives.

@Alger As suggested by Dr. Venkat Subramanian I have added implicit=true and optimize=true.

Besides that I have made a procedure that on input of the value of C and the range of t-values computes the solution.

This is convenient when different C's are considered.

Computing times are much improved. A "singularity" is still announced for C = 0.000045 and range 0..100, but this time it doesn't happen till t = 81.598332, which again makes me distrust the existence of a singularity.

MaplePrimes20-04-30_odesystemIII.mw

Using method=ck45 (i.e. Cash-Karp) instead of the default rkf45 we don't get any complaints for C=0.000045 and range = 0..100:
With the procedure res from the uploaded worksheet do:
 

VARS:=[[t,ids(t)],[t,ipdr(t)],[t,ipqr(t)],[t,iqs(t)],[t,vds(t)],[t,vqs(t)]];
dsol_45:=res(0.000045,0..100,method=ck45);
odeplot(dsol_45,VARS,0..100,numpoints=10000);
odeplot(dsol_45,VARS[-2..-1],99.9..100,numpoints=1000);

The plots:

#################################
Added May 3.

Here is a slightly revised procedure res (called res2) which doesn't require the range argument, but allows it. The intention is to find the approximate period of a potentially periodic solution.
 

#Allowing for leaving out the range simply by letting it not be a required parameter. 
#It can be added in the form range= t0..t1 in the call to res2 if desired.
res2:= proc(c::realcons) global C;
   dsolve(eval({eq1, eq2, eq3, eq4, ep5, ep6, ics},C=c), numeric, 
          maxfun = 10^8,compile=true,implicit=true,optimize=true,_rest)
end proc;
#Here the intention is to find the approximate period of a potentially periodic solution. We use #zeros of vds.
#The distance between two consequtive zeros times 2 should be the period in this case.
#We only take note of zeros beyond t=99.8 and up to 100.
dsol_45E2:=res2(0.000045,method=ck45,events=[[[vds(t),t>99.8],none]]);
dsol_45E2(99.805); # This first call takes a while
T1:=op(dsol_45E2(eventfired=[1])); # The event fired at this time
dsol_45E2(99.815); # This is fast because it starts from the result above.
T2:=op(dsol_45E2(eventfired=[1])); # The next zero
#Now automating this using the results T1 and T2:
A:=Array(datatype=float):
A(1):=T1:
A(2):=T2: 
d:=T2-T1: eps:=1e-4: # eps = tolerance
for i from 3 while T2<=100 do
  dsol_45E2(T2+d+eps);
  T1:=T2;
  T2:=op(dsol_45E2(eventfired=[1]));
  A(i):=T2;
  d:=T2-T1
end do:
  
S:=[seq(A[i+1]-A[i],i=1..numelems(A)-1)]; # The consequtive differences between zeros
listplot(S);
wr; # The period of sin(wr*t) is 2*Pi/wr
2*Pi/wr;

Why not use odeplot instead of plot in the last plots too?

Here I take the code from acer's answer:
 

CodeTools:-Usage(odeplot(dsol1,[t,(cos(wr*t)*vds(t) - sin(wr*t)*vqs(t))*sqrt(3)], 0 .. 40, labels=[t,'Uas'], color = ["Red"], style = [line],numpoints = 7500));
###
CodeTools:-Usage(odeplot(dsol1,[t,(cos(wr*t)*ids(t) - sin(wr*t)*iqs(t))], 0 .. 40, labels=[t,'Ias'], color = ["Red"], style = [line],numpoints = 7500));

Using solve:
 

restart;
eqs:=(-(2*I)*c1*eta+c1*sqrt(-e0^2-4*eta^2))/e0 = c2, ((2*I)*c2*eta+c2*sqrt(-e0^2-4*eta^2))/e0 = -c1;
solve({eqs},{c1,c2});

We see that c2 can be anything (c2=c2) and then if c2 is chosen, c1 is determined.

Go to ?pdsolve,numeric and you will find that the help page for pdsolve/numeric in the beginning in the section Parameters says:

PDEsys-single or set or list of time-dependent partial differential equations in two independent variables

Thus even if your equation had no integral pdsolve/numeric couldn't handle it.

If I understand you correctly, what you could use is the uneval modifier in your procedure.

Here is the code for simplify2 with a few changes:

restart;

simplify2 := proc (Term1::uneval)
  # Simplify a term. Additionally check if the simplification is advantageous.
  # In some cases the term gets longer after the simplification. Then discard.
  # Give the local procedure variable a long name. They must not be identical in Term.
  local simplify_tmpvar_c1, simplify_tmpvar_c1sum, simplify_tmpvar_c2, \
        simplify_tmpvar_c2sum, simplify_tmpvar_nms, simplify_tmpvar_k, Term2,Term: #Term local now
  # Assume all contents of the Term to be local variables to this procedure
  # Otherwise, variables like "l1" can be overwritten by occurences in the calling worksheet.
  Term:=eval(Term1,2); ############# change
  simplify_tmpvar_nms:=convert(indets(Term,name),list): # get list of symbols
  for simplify_tmpvar_k from 1 to LinearAlgebra:-ColumnDimension(simplify_tmpvar_nms) do
    if simplify_tmpvar_nms[simplify_tmpvar_k] = Pi or simplify_tmpvar_nms[simplify_tmpvar_k] = 0 then
      next: # 0 and Pi can not be variables
    end if:
    eval(sprintf("local %s;", String(simplify_tmpvar_nms[simplify_tmpvar_k]))); # assume local
   end do:

  # Get computational effort before simplification
  simplify_tmpvar_c1:=add(codegen:-cost~(Term)): 
  simplify_tmpvar_c1sum := diff(simplify_tmpvar_c1,functions)+ \
                           diff(simplify_tmpvar_c1,additions)+ \
                           diff(simplify_tmpvar_c1,multiplications)+\
                           diff(simplify_tmpvar_c1,divisions):
  # Perform simplification. Attention: tries to use global variables to simplify the expression (see above)
  Term2 := simplify(Term):
  # Get effort after simplification
  simplify_tmpvar_c2:=add(cost~(Term2)): 
  simplify_tmpvar_c2sum := diff(simplify_tmpvar_c2,functions)+\
                           diff(simplify_tmpvar_c2,additions)+\
                           diff(simplify_tmpvar_c2,multiplications)+\
                           diff(simplify_tmpvar_c2,divisions):
  if simplify_tmpvar_c2sum > simplify_tmpvar_c1sum then
    # The simplification was not successful. Take old version.
    Term2 := Term:
  end if:
  return Term2:
end proc:
#with(LinearAlgebra): with(ArrayTools): with(codegen): with(CodeGeneration): with(StringTools):

T := m2*(l1+l2)*(qD2)^2; # dummy-expression for kinetic energy
l1 := length(T);
T;
T_simpl := simplify2(T);

T_simpl;
l2 := length(T_simpl);

res:=[solve](1/2*m*v^2 = 1/2*F*x, v); 
res1:=simplify(res) assuming positive;
combine(res1,radical) assuming positive;

Although I do think that the exact solver in this piecewise case should do this problem correctly, it can be done manually as I have here.
It should be noted that due to a bug in Maple 2020 this only works in previous versions (well at least in Maple 2019.2 and 2018.2, and also in the one you are using Maple 2015):

restart;
with(plots):
Digits:=20:
perturbation := .9*(diff(epsilon(t), t, t)) = (-28.67085587*piecewise(0.2000000000e-1*t < 0.8e-1, 0., 0.2000000000e-1*t < .12, -0.8e-1+0.2000000000e-1*t, 0.2000000000e-1*t < .14, .16-0.2000000000e-1*t, 0.2e-1)-0.7645561571e-1*t-1.548363347)*(diff(epsilon(t), t))-0.7596363e-1*piecewise(0.2000000000e-1*t < 0.8e-1, 0., 0.2000000000e-1*t < .12, -0.8e-1+0.2000000000e-1*t, 0.2000000000e-1*t < .14, .16-0.2000000000e-1*t, 0.2e-1)-0.202569683e-3*t+0.26997403e-1-10.10*epsilon(t)
P:=simplify(perturbation,piecewise);
ivp:={P,epsilon(t0)=e0, D(epsilon)(t0)=e0_1}; # For general use
ivp0:=eval(ivp,{t0=0,e0=0,e0_1=0}) assuming t<4;
EPS1:=rhs(dsolve(ivp0));
e4,e4_1:=op((evalf@eval)~([EPS1,diff(EPS1,t)],t=4));
e4,e4_1:=op(simplify(fnormal([e4,e4_1]),zero));
ivp4:=eval(ivp,{t0=4,e0=e4,e0_1=e4_1}) assuming t>=4,t<6;
EPS2:=rhs(dsolve(ivp4));
plot(piecewise(t<4,EPS1,EPS2),t=0..6);
e6,e6_1:=op((evalf@eval)~([EPS2,diff(EPS2,t)],t=6));
e6,e6_1:=op(simplify(fnormal([e6,e6_1]),zero));
ivp6:=eval(ivp,{t0=6,e0=e6,e0_1=e6_1}) assuming t>=6,t<7;
EPS3:=rhs(dsolve(ivp6));
plot(piecewise(t<4,EPS1,t<6,EPS2,EPS3),t=0..7);
e7,e7_1:=op((evalf@eval)~([EPS3,diff(EPS3,t)],t=7)); # No small imaginary parts
ivp7:=eval(ivp,{t0=7,e0=e7,e0_1=e7_1}) assuming t>=7;
EPS4:=rhs(dsolve(ivp7));
EPS:=piecewise(t<4,EPS1,t<6,EPS2,t<7,EPS3,EPS4): # The solution
plot(EPS,t=0..10,color=blue); pexact:=%:
### Now numerically:
numsol := dsolve({perturbation, epsilon(0)=0, D(epsilon)(0)=0}, numeric,abserr=1e-20,relerr=1e-18,maxfun=0):

odeplot(numsol, [t, epsilon(t)], t=0..10, gridlines=true, color=red); pnumeric:=%:
display(pexact,pnumeric); # No difference discernible
odeplot(numsol,[t, epsilon(t)-EPS],0..10); # The difference

Note:  The problem in Maple 2020 shows up here:

EPS3:=rhs(dsolve(ivp6));
Error, (in dsolve) invalid left hand side of assignment

frontend does a good job here. Using the example by Tom Leslie (resembling yours):
 

restart;
expr:=(diff(P(t),t)+sin(Pi*x/L)*a[1](t))/(L^4*k*A*G)+sin(Pi*x/L)*diff(a[1](t),t$4)+sin(Pi*x/L)*diff(P(t),t);
####
frontend(coeff,[expr,a[1](t)]);
frontend(coeff,[expr,diff(P(t),t)]);
frontend(coeff,[expr,sin(Pi*x/L)]);