More precisely, the polygon is coloured as if it were star-shaped with respect to the first vertex.  That does not apply to the Standard interface in Maple 11, which does fill non-convex polygons correctly.  The only work-arounds I know of would be

1) take an initial vertex with respect to which the polygon is star-shaped: in your example, you could try

PLOT(POLYGONS([[0,0.4],[0,1],[1,1],[0.2,0.4],[1,0],[0,0]]),
   COLOR(RGB,1,0,0));

or 

2) break up your polygon into pieces (e.g. triangles).

The problem is basically this:

subs(x=i, F)  returns a version of F with i substituted for x, but it does not change the value of F. 

Your procedure did

subs(x = i, {delta, r, beta, alpha, m});

but then didn't do anything with the result.  The next command used delta, r, beta, alpha, m, not the values that results from substituting x=i into those expressions.

What you might try is something like this:

for i in x0 do  
  L:= [op(L), 
     OptEff(P10,P20,op(subs(x=i, [delta, beta, r, alpha, m])))]
end do;

 I assume you mean that if x[i] and x[i+1] are both variables, x[i] <> x[i+1], and similarly if x[i] and x[i+10] are both variables.

You can use integer linear programming to find a solution.

bounds:= map(t -> (t <= 9), vars) union map(t -> (t >= 1), vars);

inds:= map(op,vars);

p1:=select(t -> member(t+1,inds),inds);
p10:= select(t -> member(t+10,inds),inds);
rules:= {seq(x[i]-x[i+1]>=9*t[i]-8, i=p1),           
         seq(x[i+1]-x[i]>=1-9*t[i],i=p1),
         seq(x[i]-x[i+10]>=9*s[i]-8,i=p10),
         seq(x[i+10]-x[i]>=1-9*s[i],i=p10)};

S:=Optimization[Maximize](0,alleqs union bounds union rules,
        integervariables=vars, 
        binaryvariables={seq(t[i],i=p1),seq(s[i],i=p10)});

sort(select(has,S[2], x), (x,y)->evalb(op(lhs(x))<op(lhs(y))));

[x[22] = 7, x[23] = 9, x[26] = 8, x[27] = 7, x[28] = 9, x[32] = 9, x[33] = 8, x[35] = 8, x[36] = 9, x[37] = 5, x[38] = 7, x[42] = 7, x[43] = 9, x[44] = 6, x[45] = 9, x[46] = 4, x[53] = 4, x[54] = 3, x[56] = 5, x[57] = 3, x[64] = 1, x[65] = 6, x[66] = 1, x[67] = 6, x[68] = 2, x[72] = 9, x[73] = 7, x[74] = 4, x[75] = 1, x[77] = 2, x[78] = 3, x[82] = 2, x[83] = 3, x[84] = 1, x[87] = 1, x[88] = 2]

Warning: this can take a long time.

Is this what you mean?

with(plots):
data1:= [[0.125, -4.72557], [0.08333, -4.53944], 
  [0.0625, -4.48229], [0.05, -4.42733]];
data2:= [[0.125, -1.7312], [0.08333, -1.96658],   
  [0.0625, -2.1274], [0.05, -2.187]];
f1:= N -> -8.55704612297018*sqrt(1-0.73962210555774
      *cos(Pi/(N+1)));
f2:= N -> -8.55704612297018*sqrt(1-0.73962210555774
      *cos(Pi/(N+1)))+10.112*(1/N)+1.7295;
display([pointplot(data1),pointplot(data2),
   plot([[1/N,f1(N),N=8..20],[1/N,f2(N),N=8..20]])]);

It's impossible to tell from your description what's going wrong.  Why do you say you can't do it?  My guess is that your "intercepts" has been implicitly declared local, and you're using it as if it's global.  But that's just a guess.  The mystery could be cleared up if we could see the actual procedure definition.

I think we'll have to know the value of the variable n in order to solve this one.  It's unlikely that the problem has anything to do with the arrow command, more likely that it's from trying to evaluate n at t = 0.  Perhaps upload your worksheet?

The area between curves y = f(x) and y = g(x), for x from a to b, is

int(abs(f(x)-g(x)), x=a..b);

Whether Maple can actually find a correct exact expression for this integral is another question.  Sometimes it can and sometimes it can't.  Numerical integration  by

evalf(Int(abs(f(x)-g(x)), x = a .. b));

should work, though.

What's with the html tags in this one?

There seem to be several issues here.  One is that output from doubly-nested loops is ordinarily not printed.  For an explanation of that, see the help page ?printlevel.  But that affects only what you see, not the actual results of your code.  The reason your

for i from 1 to 3 while op(1,M[i,2])=1 do
 vars:=vars union {op(2,M[i,2])}
end do;

doesn't work is that the loop stops immediately (since op(1,M[1,2]) is not 1).

I think what you want is

for i from 1 to 3 do
  for j from 1 to 3 do
     if op(1,M[i,j])=1 then
      vars:=vars union {op(2,M[i,j])}
    end if
  end do
end do;

which ends with vars = {x[12], x[20], x[24], x[40]}.  There are simpler ways to get that,

e.g.

indets(M, name);

You can't get a polynomial for this sequence because there isn't one: the formula isn't a polynomial.  Try a formula of the form a_n = A*(-1)^n + B*n + C.

LinearMultivariateSystem should also be able to help with that.  Try this:

 S:=SolveTools[Inequality][LinearMultivariateSystem](
 {x+y+z<=1,x-y+z<=1,-x+y+z<=1,-x-y+z<=1,x>=-2,y>=-2,z>=-2},
 [x,y,z]):

makerange:= proc(S::set(relation),v::name)
   local t,a,b;
   a:= -infinity; b:= infinity;
   for t in S do
     if type(t,`=`) then 
       a:= rhs(t); b:= a
     elif lhs(t)=v then b:= rhs(t)
     else a:= lhs(t)
     end if
   end do;
   [a, b]
end proc:

makebox:= proc(L)
  plot3d(
op(subs(x=L[1,1],[[x,y,z],y=L[2,1]..L[2,2],z=L[3,1]..L[3,2]]))),
  plot3d(
op(subs(x=L[1,2],[[x,y,z],y=L[2,1]..L[2,2],z=L[3,1]..L[3,2]]))),
  plot3d(
op(subs(y=L[2,1],[[x,y,z],x=L[1,1]..L[1,2],z=L[3,1]..L[3,2]]))),
  plot3d(
op(subs(y=L[2,2],[[x,y,z],x=L[1,1]..L[1,2],z=L[3,1]..L[3,2]]))),
  plot3d(
op(subs(z=L[3,1],[[x,y,z],x=L[1,1]..L[1,2],y=L[2,1]..L[2,2]]))),
  plot3d(
op(subs(z=L[3,2],[[x,y,z],x=L[1,1]..L[1,2],y=L[2,1]..L[2,2]])));
end proc:

plots[display](map(t -> makebox(zip(makerange,t,[x,y,z])),S),
  axes=box, labels=[x,y,z]);

Note: this isn't likely to work if the feasible region is unbounded.

As of Maple 10, solve allows its arguments to be lists instead of sets, and then (I think) always returns the results with the variables in the given order:

> solve([x^2*y^2-1,x^2+y^2=1],[x,y]);

[[x = -RootOf(1+_Z^4-_Z^2)^3+RootOf(1+_Z^4-_Z^2), y = RootOf(1+_Z^4-_Z^2)], [x = RootOf(1+_Z^4-_Z^2)^3-RootOf(1+_Z^4-_Z^2), y = RootOf(1+_Z^4-_Z^2)]]

It would be nice if fsolve could be given this same capability.

I can't reproduce your problem.  One difficulty is that solve doesn't always put its results in the same order, but I think I managed to get around that.  The 0.898919459 in (28) is certainly not CSBt, it is FB.  I don't think it's a bug in Maple, more likely an editing mistake on your part: it's very easy to do some editing on Maple input and then forget to press Enter, leaving the document in a situation where some input doesn't match what appears to be its output.

That dodge of using "satisfies" would only be necessary if you needed a type check in the procedure declaration.  If you're going to do it in the body of the procedure, you can simply say

> smth := proc(T)
  if is(T,positive) then 1 else 0 end if
  end proc

As a type, "positive" is a subtype of "extended_numeric", so only integers, fractions, floats and +infinity will do.  Pi (not pi) does have the property "positive", but not the type.  If you want your procedure to accept symbolic constants, don't use "positive" in the type-checking; perhaps you could use

> ispositive := rcurry(is, positive):
  smth := proc(T::satisfies(ispositive))
...

I'm not a Matlab user, but it looks like this is supposed to solve the generalized eigenvalue problem

A x = lambda B x

using the Cholesky factorization of B.  Except for the Cholesky part, Maple's analogue would be

> LinearAlgebra[Eigenvectors](A,B);