Your worksheet is presented as a horrible mess.  I did not not have the courage to wade through it completely.  It appears that it is suffering from the 2Dmath-hormone-imbalance-syndrome.  I suggest that you consider switching to 1D input, as described in Configureing Maple, although this may be too late for you since you need to start from a fresh worksheet—I don't know how to rescue an existing worksheet.

Anyway, here is the substitution that you are looking for:

subs(
    x__s_dot(t)=diff(x__s(t),t),
    y__s_dot(t)=diff(y__s(t),t),
    ...
    u__2_dot(t)=diff(u__2(t),t),
  xEQNS):

I have shown only three of the fifteen substitutions here.  You may add the missing twelve yourself.

As we see in mehdi jafari's response, the exact solution is too complicated to be useful.  A numeric solition is more useful in this case, but for that you need to specify initial conditions, such as:

ic := A1(0)=1, D(A1)(0)=2, A2(0)=3, D(A2)(0)=4, A3(0)=5, D(A3)(0)=6;
sol := dsolve({eq1, eq2, eq3, ic}, numeric);
plots:-odeplot(sol, [[t, A1(t)],[t, A2(t)],[t, A3(t)]], t=0..6, color=[red, green, blue]);
 

By using textplot(), you may place a text anywhere you want.  For example:

restart;
with(plots):
myfont := [Times,24]:
p1 := plot(x^2, x=-1..1);
p2 := textplot([-0.2, 0.8, typeset(y=x^2)], align=left, font=myfont);
display([p1,p2], labelfont=myfont);

I can't tell what you did to get to what you have shown.  Copy and paste the following into Maple—it should do what you want:
PDE := diff(u(x,t),t) = diff(u(x,t),x,x) + 2*u(x,t)^2*(1 - u(x,t));

It is not difficult to show that the general solution of your differential equation is

where A and B are arbitrary constants.  For continuity at x=1 you need A=0, but then that leaves you with only B to play with, which means that you cannot specify two boundary conditions idependently.

Try these:

odeplot(sol, [t, Re(z(t))], t = 0 .. 3);
odeplot(sol, [t, Im(z(t))], t = 0 .. 3);
odeplot(sol, [Re(z(t)), Im(z(t))], t = 0 .. 3);

The proper way of integrating a function over a triangle is through the triangle's barycentric coordinates.  The worksheet integration-over-triangle.mw sets up a general framework for integrating any function f(x,y) over an arbitrary triangle.  At the bottom of the worksheet the general framework is applied to your specific problem.

Note added on August 22, 2017:

Robert Lopez has pointed out to me that the worksheet's computation sometimes gets the sign of the integral wrong.  This is because the worksheet assumes implicitly that the triangle's vertices are enumerated in the counter-clockwise order.  To get the correct result regardless of the ordering of the vertices, replace A with abs(A) in the expression for T_integral.

The expression which you wish to expand is a fraction whose numerator and denominator vanish with high order at [u,v]=[0,0].  A good strategy in this case is to expand the numerator and denominator separately, then combine the results, like this:

restart;
expr := ... your messy expression here ...;
Num := numer(expr):
Den := denom(expr):
N := 19;    # order of expansion
Num_N := mtaylor(Num, [u,v], N);
Den_N := mtaylor(Den, [u,v], N);
mtaylor(Num_N/Den_N, [u,v], N):  collect(%, [u,v], distributed);

The result is:

Increase N to get more terms.

Here is a self-contained worksheet for your experimentation: worksheet.mw

In your for-loop you have:

The values of all variables on the right-hand side are specified before you enter the for-loop.  Why do you expect them to change as you itererate through the loop?

Since the right-hand side is fixed, all x₁₁ have the same value, and that's why your point does not move.  You need to re-think the strategy.

An equation of the type L(n) = a*L(n-1) + b*L(n-2) is called a recurrence relationship of order two or a difference equation of order two.  To be overly pedantic, it would be called a linear homogeneous difference equation of order two with constant coefficients.

Maple's rsolve() can solve such equations explicitly.  For example:

eq := L(n) = a*L(n-1) + b*L(n-2);
sol := rsolve({eq, L(1)=p, L(2)=q}, L(n));

In your worksheet you have a=1, b=10, p=6, q=2.  We obtain the corresponding solution through:
mysol := eval(sol, {a=1, b=10, p=6, q=2});

Here we produce the first 10 terms of the sequence:
simplify([seq(mysol, n=1..10)]);
    [6, 2, 62, 82, 702, 1522, 8542, 23762, 109182, 346802]

I haven't used the FGb library myself. What I am saying here is based only on my reading of web page that you have indicated.

1. In your reply to acer, the 12627 indicates that you have downloaded an old version of the library.  Delete the dowloaded file and the files that you have extracted from it.
2. Download the newest version of the library, that is, version 1.68, which is in the first row of the table in the download page.  That will be a file named FGb-1.68.macosx.tar.gz.
3. Unpack FGb-1.68.macosx.tar.gz in the same way that you unpacked the older version.  The files will go in
        /Users/jinhuili/14539/

4. Remove the lines "mv <12627> ..." from your .mapleinit. Those are Unix commands; they don't belong to .mapleinit.

5. In your .mapleinit, change the line
       libname:= “Macintosh HD/Users/jinhuili/12627/FGblib”/FGblib, libname:
   to
       libname:= "/Users/jinhuili/14539/FGblib", libname:
   Note: I am not fully certain of the correctness of this last line because I have never used a Mac, but I am hoping that it is correct based on my experience with Linux.

6. You should be able to use the library now.  I assume that you have already installed the gmp library which is needed by FGb according to FGb's web page.

Look up "What's New" under the Help menu.

@sarra You need to name the equations as you create them so that you may refer to them later.  In MatrixAgenerating-ver2.mw I have edited your code to do that.  At the very bottom we see that you have 45 equations in 49 unknowns.  You need to supply four more equations.

restart;
with(CurveFitting):
x := [1, 2, 3, 4];
y := x^~2;
tmp := Spline(x, y, s, degree=1);
piecewise(s < x[1], y[1],
          op(tmp)[1..-2],
          s < x[-1], op(tmp)[-1],
          y[-1]);

Explanation: A[-1] is the last element of the object A.  A[-2] is the element before last.

If I understand your first question, you want an animation of the object as it rotates about the z axis.  Here is how.  I don't understand your second question.

restart;
with(PolyhedralSets):
with(plots):
with(plottools):
P1 := PolyhedralSet({
  -x-20 <= 0,
  -x-20 <= 20-y,
  -x-20 <= -(1/2)*y,
  -x-20 <= -z,
  -x-20 <= -(x+y)*(1/2),
  -x-20 <= -(y+z)*(1/2)
}):
p := Plot(P1, view = [-40 .. 40, -40 .. 40, -40 .. 40],
    color = "Orange", thickness = 2);
nframes := 60:
frames := seq(rotate(p, 0, 0, 2*Pi/nframes*i), i=0..nframes-1):
display([frames], insequence, axes=normal,
    lightmodel=light1, orientation=[70, 60, 0]);