I had a brief look at Douglas Meade's Shoot package, to which you gave a link.

It appears that in the procedure process_required_args, which is local to the module Shoot, isolate is being used like this in line 7:

ode := map(isolate,ode,diff);

Here in your case ode:=convert(ODE,list);

The problem is that the success of this isolate command depends on each differential equation not having more than one function differentiated.
Notice that your first equation does not satisfy that requirement.
The solution is to solve for the derivatives:
sys:=op(solve(ODE,diff~([F,H,f,g,u,v](eta),eta))); #Doesn't work in Maple 12
#In Maple 12:
sys:=op(solve(ODE,map(diff,[F,H,f,g,u,v](eta),eta)));
Two other problems:
1. Don't use gamma since that is Euler's constant in Maple. Change it to gm or similar.
2. FNS has H(), it obviously should be H(eta).

Then do
shoot(sys, IC, BC, FNS, [alpha = 0, gm = 0, z = -.2, Q = 0]);

But you don't get success with the guess you have.
However, following the suggestion in my comment "Shooting using dsolve/parameters" you could shoot from eta=L also when using Douglas Meade's 'shoot'. It is quite like I wrote in the comment.
Use
ICSL:={F(L) = 0, H(L) = 2, f(L) = fL, g(L) =-fL, u(L) = 0,v(L)=vL};
and then
res:=shoot(sys, ICSL, {f(0)=0,u(0)=1}, FNS, [fL=0.4,vL=0]);

Try this assuming that your actual boundary conditions are:
{F(L) = 0, H(L) = n, g(L) = -f(L), u(L) = 0,f(0) = 0,u(0)=1}

restart;
n := 2; M := 3; L := 6; B := 0.2e-1;
ODE := {g(eta)*(diff(g(eta), eta))+B*(f(eta)+g(eta)) = 0, g(eta)*(diff(F(eta), eta))+F(eta)^2+B*(F(eta)-u(eta)) = 0, g(eta)*(diff(H(eta), eta))+H(eta)*(diff(g(eta), eta))+F(eta)*H(eta) = 0, diff(v(eta), eta)+f(eta)*v(eta)-u(eta)^2+B*H(eta)*(F(eta)-u(eta))-M*u(eta) = 0, diff(f(eta), eta) = u(eta), diff(u(eta), eta) = v(eta)};
BC := {F(L) = 0, H(L) = n, g(L) = -f(L), u(L) = 0,f(0) = 0,u(0)=1};
sol:=dsolve((ODE union BC),numeric,approxsoln=[F(eta)=0,H(eta)=2,f(eta)=1,g(eta)=-.2,u(eta)=0,v(eta)=0]);
plots:-odeplot(sol,[[eta,10*F(eta)],[eta,H(eta)],[eta,f(eta)],[eta,g(eta)],[eta,u(eta)],[eta,v(eta)]],0..L);

Or plotting the solution as an array:

A:=[[eta,F(eta)],[eta,H(eta)],[eta,f(eta)],[eta,g(eta)],[eta,u(eta)],[eta,v(eta)]];
plots:-display(Array([[seq(plots:-odeplot(sol,A[i],0..L),i=1..3)],[seq(plots:-odeplot(sol,A[i],0..L),i=4..6)]]));

Having seen the output from Maple and realized that it could be written as a linear combination of cosh terms, I did the following, which gives a much neater solution:

restart;
eq:=a11*diff(phi(x),x,x,x,x)+a22*diff(phi(x),x,x)+a33*phi(x)=0;
bcs:=phi(a)=sigma1 , phi(-a)=sigma1 , D(phi)(a)=0 , D(phi)(-a)=0;
#The characteristic polynomial
pol:=a11*r^4+a22*r^2+a33; #quadratic in r^2 (!)
#The actual parameters:
params:={a11=2.731e-10,a22=-1.651e-9,a33=3.09027e-10,a=35.714,sigma1=200e6};
fsolve(eval(pol,params),r);
#We see that we are in the 4 real roots case (2 positive and 2 negative).
#Thus the symmetry of the bcs means that there will be a solution on the form
sol:=A*cosh(r1*x)+B*cosh(r2*x);
# where r1 and r2 are two roots of pol with different absolute value.
#Finding A and B from conditions at one end. Symmetry does the rest:
eqs:={eval(sol=sigma1,x=a),eval(diff(sol,x)=0,x=a)};
AB:=solve(eqs,{A,B});
R:=solve(pol=0,r); #The roots
R[1]+R[2]; #Gives 0, so abs-values are equal
R[3]+R[4]; #Same
sol2:=subs(AB,r1=R[1],r2=R[3],sol); #R[1] and R[3] have different abs-value.
odetest(phi(x)=sol2,[eq,bcs]); #Checking, 5 zeros: OK
u:=evalf(eval(sol2,params));
plot(u,x=-35.714..35.714);

I solve without using the actual values first. I make the same corrections as Kitonum does in his answer.

restart;
eq:=a11*diff(phi(x),x,x,x,x)+a22*diff(phi(x),x,x)+a33*phi(x)=0;
bcs:=phi(a)=sigma1 , phi(-a)=sigma1 , D(phi)(a)=0 , D(phi)(-a)=0;
params:={a11=2.731e-10,a22=-1.651e-9,a33=3.09027e-10,a=35.714,sigma1=200e6};

sol:=dsolve({eq,bcs},phi(x)); #You need the second argument in this non concrete case
length(sol); #A measure of the size of the output: Huge do to solving a quartic.
odetest(sol,[eq,bcs]); #5 zeros means OK
u:=simplify(rhs(sol),size); #The result can be shortened
length(u);
Digits:=20: #Since rounding errors are inevitable raising Digits is necessary
u1:=evalf(eval(u,params));
#Checking if Digits may be high enough for bcs
eval(u1,x=-35.714);
eval(u1,x=35.714);
eval(diff(u1,x),x=-35.714);
eval(diff(u1,x),x=35.714);
plot(u1,x=-35.714..35.714);

Use 1D input (a.k.a. Maple input). Change preferences in Tools/Options/Display.

And while at it, why not change to Maple Worksheet in Tools/Options/Interface.

Example:

A:=Matrix([[1,2],[3,4]]); #Evidently not symmetric
B:=Matrix(A,shape=symmetric); #This is

I think you have an implied multiplication in your 2D math input:

assign (space) (L[2]);

I would just insert alpha and beta into the equations:
sys:=diff(y(t),t)=alpha0*(x(t)-x0)-y(t), diff(x(t),t)=y(t)-beta0*(x(t)-x0);
#The general solution:
res:=dsolve({sys});
#The solution satisfying x(0)=x0, y(0)=y0
sol:=dsolve({sys,x(0)=x0,y(0)=y0});
#Extracting x(t) and y(t):
X,Y:=op(subs(sol,[x(t),y(t)]));
#Plotting example:
plot(eval([X,Y],{x0=.1,y0=0.2,alpha0=1.2345,beta0=.6789}),t=0..5);

It is just as easy to plot alpha0*(X-x0) and beta0*(X-x0).

(Did you really want alpha and beta both to have x(t)-x0 as a factor?)

Shouldn't eq1 and eq2 be:
eq1 := .5*r*(u-sin(u)) = b;
eq2 := a-.5*r*(1-cos(u)) = 0;

And what about the minus sign in the square root in the integral?

Have a look at Student:-Calculus1:-RiemannSum, specifically the partition option and the method option.

Illustration:

# Illustrating the method=procedure option with a rather simple example and in an abstract setting
restart;
p:=proc(f,x,c,d) c+(d-c)/3 end proc;
Student:-Calculus1:-RiemannSum(f(x),x=a..b,method=p,partition=3);
#Illustrating the partition option using partition=list of points:
n:=5:
pt:=[a,seq(a+(b-a)/n*k+(b-a)/n^2,k=0..n-1),b]; #Forgot a first time around.
Student:-Calculus1:-RiemannSum(f(x),x=a..b,method=p,partition=pt);
#Making f, a, and b concrete:
f:=x^2:
a:=2: b:=5:
pt;
Student:-Calculus1:-RiemannSum(f,x=a..b,method=p,partition=pt);
evalf(%);
int(f,x=a..b);

Any solution to your linear ode is a particular solution. Thus you only need to test whether the solution found is indeed a solution:

odetest(m,ode); #Should return 0, and does!

If some other answer is given (also correct, like the one you cite) then the difference between the two will solve the corresponding homogeneous equation.

The second solution has the advantage that is shorter. That is all.

If you want to find just one solution the easiest way is (in this case) simply to find the general solution

dsolve(ode);

and then set the arbitrary constants _C1 and _C2 to any constant value you like , e.g. 0.

The error message says that you cannot have boundary conditions with derivatives not normal to the boundary.
In your case the problem is the condition D[2](u)(0, t) = D[2](u)(2, t), which involves the time derivatives (the number 2) at the boundaries x=0 and x = 2.

Maybe what you actually meant was  D[1,1](u)(0, t) = D[1,1](u)(2, t), which involves the second derivatives w.r.t. the first variable x.
 

Your initial conditions have z(0)=0, thus k(0) will result in a division by zero. Change z(0) to any number <> 0.

Test outside procedure:

N:=10: #N clearly has to be concrete.
F := dsolve(eval({q1,q2,q3,q4,ic},{a=5,b=7}),numeric, output=Array([seq(k,k=0..N)]));
plots:-odeplot(F,[t,x(t)],0..N);

If I understand you correctly, you could do like this starting with your definition of the lists (omitted here):

Dim := <L, M, Q, R, T, theta>;
A:=Matrix(numelems(Dim),8);
for j to 8 do A[..,j]:=evalindets(subs(map(`=`@op,l[j]),Dim),name,0) end do:
A;
NoOflists := Vector[row]([``,seq(list[i], i = 1 .. 8)]);
<NoOflists,<Dim|A>>;

#Comment: It seems also to work if e.g. l[8]:=[]; (The empty list), but not if l[8] is not assigned. 

The correct syntax for

`&PartialD;`(`&PartialD;`(W(x, y))/`&PartialD;`(x))/`&PartialD;`(x);

is

diff(W(x,y),x,x);

Similarly the correct syntax for

`&PartialD;`(`&PartialD;`(Z(x, y))/`&PartialD;`(x))/`&PartialD;`(y);

is

diff(Z(x,y),x,y);