@ssteve It is confusing that you unapply Lfinal with respect also to  x since Lfinal doesn't depend on x.
Consequently m is really just a function of 2 variables.

Here is a start. The analysis is your job.

restart;
ode:=diff(x(t),t,t)+b*diff(x(t),t)+sin(x(t))=F*cos(x(t));
res:=dsolve({eval(ode,{b=.22,F=2.7}),x(0)=1,D(x)(0)=0},numeric,output=listprocedure,maxfun=10^6);
plots:-odeplot(res,[x(t),diff(x(t),t)],0..25);
plots:-odeplot(res,[x(t),diff(x(t),t)],1000..2000,numpoints=10000);
X,X1:=op(subs(res,[x(t),diff(x(t),t)]));
M:=Matrix([seq([X(k*2*Pi),X1(k*2*Pi)],k=100..5000)]);
plot(M,style=point);
#Or use poinplot
plots:-pointplot(M);

u:=sin(w*t-theta)+sin(w*t-theta-2*Pi/3)-sin(w*t-theta+2*Pi/3);
expand(%);
combine(%);
convert(%,phaseamp,t);

If I understand your intention correctly here is an example where for each value of the loop variable a column is recorded in a matrix:

restart;
N:=5: d:=3:
M:=Matrix(d,N):
for i from 1 to N do
   M[1..d,i]:=LinearAlgebra:-RandomVector(d)
end do:
M;

restart;
with(LinearAlgebra):
J1:=2:
myGamma:=0.8: lambda12:=2: lambda21:=2: eta12:=2: eta21:=2:
R0:= myBeta/(J1*myGamma);
u2:=1-v1-v2-u1:
eq1:=-eta12*v1+eta21*v2+myBeta*u1*v1-myGamma*v1:
eq2:=-eta21*v2+eta12*v1+myBeta*u2*v2-myGamma*v2:
eq3:=-lambda12*u1+lambda21*u2-myBeta*u1*v1+myGamma*v1:
sys:=[eq1,eq2,eq3];
var:=[v1,v2,u1]:
J:=unapply(VectorCalculus:-Jacobian(sys, var),var):
res:=solve(sys,{v1,v2,u1});
nops([res]);
indets([res],RootOf);
op([1,1],%);
discrim(%,_Z);
solve(%);
#Thus if myBeta is between those values the RootOf-point has imaginary components.
#Assuming that you are interested only in -32/5 < myBeta < 48/5 we got the following 2 points:
res1:=remove(has,[res],RootOf);
#One is independent of myBeta
pt1:=[0,0,0.5];
Eigenvalues(J(op(pt1)));
#We see that that point is asymptotically stable if  myBeta < 1.6 and unstable if myBeta > 1.6.
#Now the other point:
select(has,res1,myBeta);
pt2:=subs(op(%),var);
Eigenvalues(J(op(pt2)));
#We see that the second point is asymptotically stable if  myBeta > 1.6 and unstable if    -32/5< myBeta < 1.6.
pt24:=[op(pt2),simplify(eval(u2,var=~pt2))];
#Notice that v2 = v1 and u2 = u1, so the plot will show only two graphs (and colors will be chosen accordingly):
plot(pt24,myBeta=1.6..48/5,linestyle=1,color=[red,red,blue,blue]):
plot(pt24,myBeta=1..1.6,linestyle=2,color=[red,red,blue,blue]):
plots:-display(%%,%); p2:=%:
pt14:=[op(pt1),eval(u2,var=~pt1)];
plot(pt14,myBeta=1.6..48/5,linestyle=2,color=[maroon,maroon,green,green],thickness=2):
plot(pt14,myBeta=1..1.6,linestyle=1,color=[maroon,maroon,green,green],thickness=2):
plots:-display(%%,%); p1:=%:
plots:-display(p1,p2);

On the picture I see a dot after the summation sign. If that is a multiplication sign then it doesn't belong there.

Also I see square brackets. Those cannot be used for grouping in Maple as they are used for a lot of other things.
Use ordinary parentheses.

In w a and b are the formal parameters used in the definition. The value of g contains the global variables a and b.
So change w to something like this:

w:=proc(aa,bb) eval(g,{a=aa,b=bb})+14 end proc;

The problem is that diff(theta(t),t) is seen in diff(theta(t),t,t) = diff(diff(theta(t),t),t).
So you could use frontend which will freeze expressions like that:

frontend(coeff,[rhs(eq1pp),diff(theta(t),t)]);

Although frontend is not necessary for the second derivative it works for that too:

frontend(coeff,[rhs(eq1pp),diff(theta(t),t,t)]);

Use seq instead:

p:=2*Vector[column]([seq(i,i=Nr..0,-1)])+LinearAlgebra:-ConstantVector[column](1,Nr+1);

For an example I create a list of lists randomly. I haven't bothered to make sure one is positive the other negative.

L:=RandomTools:-Generate(listlist(float(range=-1..1),7,2));
#Selecting positive elements. In your case the result will be a list of lists with each just one element.
select~(type,L,positive);
#Removing the list brackets at the same time:
L1:=(op@select)~(type,L,positive);
#In your case L1 will have as many elements as L
#Finding the minimum
min(op(L1));
#Finding the minimum and its position (if that is important)
ListTools:-FindMinimalElement(L1,position);

For the stability you want to look at the eigenvalues of the jacobian.

From your worksheet it is not clear if  the left hand sides of eq1,eq2,eq3 corresponds to v1' ,v2', u1' or some other ordering. I shall assume that that is the ordering:

J:=unapply(VectorCalculus:-Jacobian(lhs~([eq1,eq2,eq3]),[v1,v2,u1]),v1,v2,u1);
J(0,0,1/2);
LinearAlgebra:-Eigenvalues(%);

If the ordering is different then just change the order v1,v2,u1  in two places in the definition of J.

Addendum: In response to a question on how to plot the equilibrium points for a range of myBeta-values:

restart;
J:=2:
myGamma:=0.8: lambda12:=2: lambda21:=2: eta12:=2: eta21:=2:
R0:= myBeta/(J*myGamma);
u2:=1-v1-v2-u1:
eq1:=-eta12*v1+eta21*v2+myBeta*u1*v1-myGamma*v1:
eq2:=-eta21*v2+eta12*v1+myBeta*u2*v2-myGamma*v2:
eq3:=-lambda12*u1+lambda21*u2-myBeta*u1*v1+myGamma*v1:
sys:=[eq1,eq2,eq3];
var:=[v1,v2,u1]:
#Make a procedure for the equilibrium points:
Pts:=proc(beta) local SYS,sol,i,p;
   if not type(beta,numeric) then return 'procname(beta)' end if;
   SYS:=eval(sys,myBeta=beta);
   sol:=RootFinding:-Isolate(SYS, var, digits = 4);
   p:=NULL;
   for i from 1 to  nops(sol) do
      p := p,subs(sol[i],var)
   end do;
   p
end proc;
#----------------------------------
#If desired include u2:
VU:=proc(beta) local SYS,sol,i,p;
   if not type(beta,numeric) then return 'procname(beta)' end if;
   SYS:=eval(sys,myBeta=beta);
   sol:=RootFinding:-Isolate(SYS, var, digits = 4);
   p:=NULL;
   for i from 1 to  nops(sol) do
      p := p,subs(sol[i],[v1,v2,u1,u2])
   end do;
   p
end proc;
#Testing:
Pts(1.9);
VU(1.9);
#Animating the equilibrium points:
plots:-animate(plots:-pointplot3d,[[Pts(beta)],symbolsize=20,symbol=solidsphere,axes=boxed,labels=var],beta=1..4);
#Plotting v1,v2,u1,u2:
plot([seq(VU(beta)[1][i],i=1..4)],beta=1..4,color=[red,blue,green,maroon]);
plot([seq(VU(beta)[2][i],i=1..4)],beta=1..4,color=[red,blue,green,maroon]);

One problem with the integral is that its integrand contains the product of the second derivatives of y[A] and alpha[1].
Thus adding eqn[10] makes the whole system much more complicated.
The integrand does have a special structure, which could be used.
It is h(t)^2*diff(h(t),t)^2, where h is given in heqn below.
Including heqn with elim[2] (without eqn[10]) you could do like this, where ics9 leaves out G(0)=0:

heqn:=h(t)=diff(y[A](t), t)-130*sin(alpha[1](t))*(diff(alpha[1](t), t));
sol9 := dsolve(elim[2] union {ics9} union {heqn}, numeric, maxfun = 10^6,output=listprocedure);
H:=subs(sol9,h(t));
H(.1); #Test
plots:-odeplot(sol9,[t,h(t)],0..2);
#Now the idea is to use fdiff and numerical integration:
fdiff(H,[1],[.1]); #Test
g:=proc(t) H(t)^2*fdiff(H,[1],[t])^2 end proc; #The integrand
g(.1); #Test
res:=evalf(Int(g,0..2,epsilon=1e-3));
# (I left out the constant 68.89).

Would this do it:

S,R:=selectremove(has,eq2,phi[i](y)*diff(w(y),y));
Parts(S,phi[i](y));
eq3:=R+%;

There is a problem in isolating diff(T(y),y,y) in sys[2] since a quadratic appears, thus there is ambiguity.

Try

sysT:=solve(sys[2],{diff(T(y),y,y)});

You get 2 solutions. Choose one of them say the first:

p1:=dsolve({sys[1],op(sysT[1]),sys[3..4],f(0)=tw,D(f)(0)=1,T(0)=1,S(0)=1,G(0)=1,D(f)(6)=0,T(6)=0,S(6)=0,G(6)=0},fcns,type=numeric,method=bvp[midrich],abserr=1e-10);

Now you get an error message saying that there are too few boundary conditions: expected 10, got 9.
So you have to deal with that.

restart;
A := Matrix( [ [a*s, a*t, a*r, a*l],[b*s, b*t, b*r, b*l],
                    [c*s, c*t, c*r, c*l], [d*s, d*t, d*r, d*l] ] );

X:=Matrix(4,symbol=x); #No functional notation
eq1 := a*d = b*c, s*l = t*r;
U := convert( Equate(A, X), set );
res := SolveTools:-PolynomialSystem(U union {eq1},indets(U union {eq1}), maxsols=1);

I noticed that the syntax has changed since Maple12. This syntax works in Maple 16 and so does the shorter:

res := SolveTools:-PolynomialSystem(U union {eq1}, maxsols=1);