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);

I looked at wikipedia:
http://en.wikipedia.org/wiki/Successive_over-relaxation

and wrote code following the formula from that link:

Using your A and b:
A:=Matrix(evalf(A),datatype=float);
b:=Vector(evalf(b),datatype=float);
n:=8:
omega:=1.136:
x:=Vector(n,0.,datatype=float): #Start value for x
Nmax:=30000: #Number of iterations
to Nmax do
  for i from 1 to n do
    x[i]:=omega/A[i,i]*(b[i]-add(A[i,j]*x[j],j=1..i-1)-add(A[i,j]*x[j],j=i+1..n))-(omega-1)*x[i];
  end do
end do:
x;

It seems that omega = 1.136 is OK, but that you cannot get much higher without divergence.