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

You have posted a screenshot of your worksheet.  That's not a good idea.  You should really post the worksheet itself to allow people to make sensible comments about it.

If you post a question in the future, note the big fat green up-arrow in the compose area.  Click on that arrow to upload your worksheet.

Anyway, here is an example of how to solve a system of linear equations in Maple.  Suppose we want to solve the system M*x = N with some 2x2 matrix M and a 2-vector N.  We do:

restart;
with(LinearAlgebra):
M := < 1,2; 3, 4>;
N := < 5, 6>;
LinearSolve(M,N);

This will produce the solution <-4, 9/2> of the system.

I haven't used Maple's pdsolve() so I cannot answer your question.  I hope that others will.  I just want to point out that the solution that you have obtained is incorrect.  For one thing, the solution should not oscillate, and certainly it should not go below zero because that would violate the maximum principle.

It is likely that pdsolve's help page discusses ways to increase the solver's accuracy.  Have a look.

Add the discont=true option, as in:

plot(tan(x), x=-2*Pi..2*Pi, discont=true);

The function MatrixExponential(A,t) in the LinearAlgebra package yields the fundamental matrix corresponding to any square matrix A.

The fundamental matrix plays an important conceptual role in the theory of linear systems of differential equations but it's not particularly useful as a computational tool. Do MatrixExponential(A,t) with your matrix A to see what I mean.
 

Aside: You have typed "eigenvectors(A)".  You probably know that Maple is case-sensitive.  The correct command is "Eigenvectors(A)" with uppercase letter E.