Try this method first. After that you get off easier by the method at the bottom.
1). Solve the Eq1 bvp problem first.
2). Notice that if f solves that problem then g(x) = -f(-x) (x in -1..0) will solve Eq1 with g(-1)=0, D(g)(-1)=1, g(0)=0,(D@D)(g)(0)=0, D(g)(0)=D(f)(0), and (D@@3)(g)(0) = (D@@3)(f)(0).
3). It follows that we can extend the f found in 1) smoothly by f(x) = -f(-x) for x in -1..0.
4). Now solve the bvp problem for Eq2 with the known option.

Eq1:=diff(f(x),x$4)+(f(x)*diff(f(x),x$3)-diff(f(x),x$1)*diff(f(x),x$2))=0;
Eq2:=diff(f2(x),x$2)+f(x)*diff(f2(x),x$1)=0;
bcsf:=f(0)=0,(D@@2)(f)(0)=0,f(1)=0,D(f)(1)=1;
resf:=dsolve({Eq1,bcsf},numeric,output = listprocedure);
F,F1,F2,F3:=op(subs(resf,[seq(diff(f(x),[x$i]),i=0..3)]));
##Illustrations of smoothness:
plot(-F(-x),x=-1..0):
plot(F,0..1):
plots:-display(%%,%);
plot(F1(-x),x=-1..0):
plot(F1,0..1):
plots:-display(%%,%);
plot(-F2(-x),x=-1..0):
plot(F2,0..1):
plots:-display(%%,%);
plot(F3(-x),x=-1..0):
plot(F3,0..1):
plots:-display(%%,%);
##The extended f:
FF:=x->piecewise(x<0,-F(-x),F(x));
resf2:=dsolve({subs(f=FF,Eq2),f2(-1)=0.5,f2(1)=1},numeric,output = listprocedure,known=FF);
plots:-odeplot(resf2,[x,f2(x)],-1..1);
##############
##Now that we know that f is going to be an odd function we can solve the problem much easier:
bcsfB:=f(-1)=0,(D)(f)(-1)=1,f(1)=0,D(f)(1)=1,f2(-1)=0.5,f2(1)=1;
resfB:=dsolve({Eq1,Eq2,bcsfB},numeric,output = listprocedure);
plots:-odeplot(resfB,[[x,f(x)],[x,f2(x)]],-1..1);

Maybe I misunderstand your problem because it seems obvious to me: You need
alpha>=-(1/2)*S^2/(K+2)
in order to avoid imaginary numbers.
Plotting should not in principle be a problem either since plot just ignores imaginary numbers:
plot3d(eval(lambda1,K=0.12345),alpha=-15..15,S=-15..15);
plot(eval(lambda1,{K=0.12345,S=2.678}),alpha=-15..15);

The following version of your first method seems to give the correct coefficients:

AdamsMoulton := proc (k::posint)
  local P,x,f,y,n,t,h,ynp1,cfs,v;
  P := CurveFitting:-PolynomialInterpolation([seq(t[n]-j*h, j = -1 .. k-2)], [seq(f[n-j], j = -1 .. k-2)], x);
  ynp1:=simplify(int(P, x = t[n] .. t[n]+h));
  cfs:=coeffs(eval(ynp1,h=1),[seq(f[n-j],j=-1..k-2)],v);
  subs([v]=~[cfs],[seq(f[n-j],j=-1..k-2)])
end proc:

AdamsMoulton(5);

Notice that 'interp' has been deprecated. Thus I used CurveFitting:-PolynomialInterpolation instead.

I found the following nice and concise reference:

http://www.math.odu.edu/~jhh/chs8.pdf

apparently written by John H. Heinbockel  from Old Dominian University in Virginia.

I wrote a Maple procedure following his algorithm.
No need to restrict to autonomous equations, so I didn't.
See the algorithm for explanations.
It seems to work fine!

TMxy:=proc(f,x0,y0,h,xn,N::posint) local x,y,ode,u,i,Tn,res,X0,Y0;
   ode:=diff(y(x),x)=f(x,y(x));
   u[1]:=f(x,y(x));
   for i from 2 to N do
      u[i]:=eval(diff(u[i-1],x),ode)
   end do;
   Tn:=unapply(eval(add(u[i]*h^(i-1)/i!,i=1..N),y(x)=y),x,y); #f depends on x and y
   res[0]:=[x0,y0];
   X0:=x0;
   Y0:=y0;
   for i from 1 to ceil((xn-x0)/h) while X0<=xn do
      Y0:=Y0+h*Tn(X0,Y0);
      X0:=X0+h;
      res[i]:=[X0,Y0]
   end do;
   [seq(res[j],j=0..i-1)]
end proc;
#Example:
Exact:=dsolve({diff(y(x),x)=x+y(x),y(0)=1});
res:=TMxy((x,y)->x+y,0,1,0.1,2,4);
p1:=plot(res,color=red):
p2:=plot(rhs(Exact),x=0..2):
plots:-display(p1,p2);

Your last point is a wrong argument.
Take the example product(exp(-1/n),n=1..infinity);
Clearly exp(-1/n) -> 1 as n -> infinity (increasingly).
However, taking logarithms we consider

-sum(1/n,n=1..infinity) = - infinity.
Thus product(exp(-1/n),n=1..infinity) = 0.

In your case, though, using cos(x) >= 1-x^2/2 (for small x) combined with ln(1-y) >= -2*y for small y>0, we get for some large enough N:

sum(ln(cos(Pi/k)),k=N..infinity) >= - sum( (Pi/k)^2, k=N..infinity) > -infinity.
Thus the product you consider is certainly not zero.

(In fact the cosine inequality holds for all x and the ln(1-y) inequality holds for 0<= y <= 0.7968 .. and y = (Pi/k)^2/2 = 0.5483 .. for k = 3. Thus you may take N=3 above).

You should not omit x in Tg, Tw, Tz.
Thus the equations should be

eqa1 := A1*(diff(Tg(x), x, x))+A2*(diff(Tg(x), x))+(A3+A4)*Tg(x)+A3*Tz(x)+A4*Tw(x) = 0;

eqa2 := B1*(diff(Tw(x), x, x))+B2*(diff(Tw(x), x))+(B3+B4)*Tw(x)+B3*Tz(x)+B4*Tg(x) = 0;

eqa3 := C1*(diff(Tz(x), x, x))+(C3+C4)*Tz(x)+C3*Tg(x)+C4*Tw(x) = 0;

#And
dsolve({eqa1,eqa2,eqa3});

#returns a large and complicated looking result.
#This version returns a more compact looking result:

dsolve({eqa1,eqa2,eqa3},{Tg(x),Tw(x),Tz(x)},method=laplace);
#Notice in this case the need for the second argument (the unknown functions).

##  RootOf has to be used since there are no formulas for the roots of polynomials of degree 5 or higher. Your polynomial is of degree 6.

You mention boundary conditions. I don't see then in your posted question.
You may want to solve numerically instead, since not much insight is gained from looking at a very complicated expression. For that you need concrete values for all constants.
See an example below.

Version 2 of your system should solve without any problem, but is still very complicated and (of course) still involves roots of a polynomial of degree 6.

#######Numerical solution example:
indets({eqa1,eqa2,eqa3},name) minus {x};
params:= %=~{seq(i,i=1..nops(%))}; #My arbitrary choice of constants.
res:=dsolve(eval({eqa1,eqa2,eqa3},params) union {Tg(0)=0,D(Tg)(0)=1,Tw(0)=1,D(Tw)(0)=0,Tz(0)=0,D(Tz)(0)=0},numeric);
plots:-odeplot(res,[[x,Tg(x)],[x,Tw(x)],[x,Tz(x)]],0..10);

RK3((x,y)->-y,0,h,1,1);
rk3:=expand(%[2,2]);
taylor(eval(rhs(res),x=h),h=0);
taylor(eval(rhs(res),x=h)-rk3,h=0);

Thus the error in a single step is O(h^4).
Since h = (b-a)/N the error in taking N steps will be N*O(h^4) = ((b-a)/h)*O(h^4) = O(h^3).

If you want a plot:

Digits:=30:
plots:-loglogplot(abs(eval(rhs(res),x=h)-rk3),h=1e-5..1);

sys:={diff(x(t),t,t)=0,diff(y(t),t,t)=0, x(0)=0.5, D(x)(0)=0.2, y(0)=0.5, D(y)(0)= sqrt(3)/5};
res:=dsolve(sys,numeric,events=[[x(t)=0,diff(x(t),t)=-diff(x(t),t)],[x(t)=1,diff(x(t),t)=-diff(x(t),t)],[y(t)=0,diff(y(t),t)=-diff(y(t),t)],[y(t)=1,diff(y(t),t)=-diff(y(t),t)]]);
plots:-odeplot(res,[x(t),y(t)],0..20,axes=boxed,frames=100);

Adding the range argument:

res:=dsolve(sys,numeric,events=[[x(t)=0,diff(x(t),t)=-diff(x(t),t)],[x(t)=1,diff(x(t),t)=-diff(x(t),t)],[y(t)=0,diff(y(t),t)=-diff(y(t),t)],[y(t)=1,diff(y(t),t)=-diff(y(t),t)]],range=0..100);
plots:-odeplot(res,[x(t),y(t)],0..100,axes=boxed,refine=1);
plots:-odeplot(res,[x(t),y(t)],0..100,axes=boxed,frames=500,numpoints=1000);

You may try statementwise as in
Matlab:-FromMatlab("y(1,:) = ya");

which prints
y(1, ..) := copy(ya);

You can look at the help page

?dsolve,numeric,ivp
ode:=diff(x(t),t)=sin(t*x(t));
res:=dsolve({ode,x(0)=1},numeric,method=classical[abmoulton]);
plots:-odeplot(res,[t,x(t)],0..10);

That doesn't give you the code, though.

Also take a look at

?Student:-NumericalAnalysis

I made my own (perfect) data in this example:

f:=x=a*y^b*z^c;
Y:=Vector(21,i->(i-1)/20.);
Z:=Vector(21,i->2+(i-1)/20.);
eval(rhs(f),{a=3,b=2,c=-3});
X:=Vector(21,i->3*Y[i]^2/Z[i]^3);
A:=Matrix([Y,Z,X],datatype=float);

Statistics:-NonlinearFit(rhs(f),A,[y,z]);

This would do the job:
restart;
plot(piecewise(x>=0,x^2,undefined),x=-2..2);

I don't think that this is a bug, but rather a weakness or limitation.

You could do this after the output with the list of severeal solutions:

select(x->op([1,2],x)>=0,%);

But notice that surely you don't have all solutions satisfying omega[n]>=0.

To illustrate that do:
expand(eqs);
subs(%[1],%[2]);
EQ:=isolate(%,cot(omega[n]));
plot([lhs,rhs](EQ),omega[n]=-10..10);

I had no problem with this:

restart;
r:=exp(x*I*phi)+x;
i := 0;
p := Array(0 .. 10);
for x from 0 by .2 to 2
  do p[i] := plots:-complexplot(r, phi = -(1/2)*Pi .. (1/2)*Pi);
  i := i+1
end do;

p[10];
p[1];
plots:-display(p[1],p[10]);
plots:-display(seq(p[i],i=0..10),insequence=true); #Animation

If you have problems, you should show us how r is defined.

Edited.
restart;
ode:=-diff(psi(x),x,x)+V(x)*psi(x)=E*psi(x);
res:=dsolve(eval(ode,V(x)=K0)); #v(x)=K0 for x<=0.
collect(convert(res,exp),exp);
#We see that k0=sqrt(E-K0). Similarly kL=sqrt(E-KL), where V(x) = KL for x >=L 
#We try with a simple potential:
Vx:=piecewise(x<=0,K0,x<L,(KL-K0)/L*x+K0,KL);
plot(eval(Vx,{L=1,K0=1,KL=2}),x=-2..4);
#The boundary value problem will be solved numerically and needs to be split in real and imaginary parts:
ode1:=eval(ode,psi(x)=u(x));
ode2:=eval(ode,psi(x)=v(x));
#The boundary conditions:
eval( psi(0)-I*D(psi)(0)/sqrt(E-K0)=2*a0,{psi=u+I*v,a0=Re(a0)+I*Im(a0)});
expand(%);
c1:=map2(remove,has,%,I);
c2:=expand((%%-c1)/I);
eval(psi(L)+I*D(psi)(L)/sqrt(E-KL)=2*aL,{psi=u+I*v,aL=Re(aL)+I*Im(aL)});
expand(%);
c3:=map2(remove,has,%,I);
c4:=expand((%%-c3)/I);
#Making the boundary conditions functions of L,K0,KL,a0,aL,E:
bcs:=unapply({c1,c2,c3,c4},L,K0,KL,a0,aL,E);

#Solving on -infinity..infinity piecewise.
#The first argument is the given potential (as a function):
#The last argument is E but needs another name, since E is present as a global variable in ode1 and ode2.
res:=proc(W,L,a0,aL,EE) local sol,u0,u1,v0,v1,k0,kL,b0,bL;
   sol:=dsolve(eval({ode1,ode2},{V=W,E=EE}) union bcs(L,W(0),W(L),a0,aL,EE),numeric,output=listprocedure);
   u0,u1,v0,v1:=op(subs(sol,[u(x),diff(u(x),x),v(x),diff(v(x),x)]));
   k0:=sqrt(EE-W(0));
   kL:=sqrt(EE-W(L));
   b0:=1/2*(u0(0)+I*v0(0)+I*(u1(0)+I*v1(0))/k0);
   bL:=1/2*(u0(L)+I*v0(L)-I*(u1(L)+I*v1(L))/kL);
   piecewise(x<0,a0*exp(I*k0*x)+b0*exp(-I*k0*x),x<=L,u0(x)+I*v0(x),aL*exp(I*kL*(L-x))+bL*exp(-I*kL*(L-x)))
end proc;
W:=unapply(eval(Vx,{K0=1,KL=2,L=1}),x); #Making the potential very concrete.
plot(W,-2..3);
sol:=res(W,1,3+2*I,5,3); #Example, input W, L, a0, aL, E.
plots:-complexplot(sol,x=-5..5); #Plot of the solution in the complex plane
plot(Re(sol),x=-5..5);
plot(Im(sol),x=-5..5);