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.

Follow these steps:

1. Sketch a graph of the function defined in your equation (6).
2. Did you do that?
   If not, then go to Step 1.
3. You DID sketch a graph of the function defined in equation (6), right?
   Now go to step 4.
4. What is the area under that graph?  You should be able to tell that the area is zero by just looking at it.  No computation is necessary.
5. Done!

Aside: Where you say Kronecker delta you should say Dirac delta.  Kronecker delta is something quite different.  What you have is neither.

In your worksheet you have:

x := seq(x =-1..1, dx);

That doesn't make much sense.  Change that to:

x := Array(0..nx, i-> 2/nx*i - 1, datatype=float[8]):

Afterward, when you want to refer to the entries of x, use square brackets, as in x[i].

There is no solution in terms of elementary functions, but you can obtain a solution numerically if you provide intitial conditions, as in:

de := -1.699678499+diff(z(t),t,t) - 3.813322790*z(t) + 0.6966019659*z(t)^2 + 0.01548004368/z(t)^2 = 0;
ic := z(0)=1, D(z)(0)=0;
dsol := dsolve({de,ic}, numeric);
plots:-odeplot(dsol, t=0..8);

which produces the following plot:

Note that in specifying your differential equation, every z should be entered as z(t).

Have a look at this worksheet.  It does some calculations with your problem but I don't quite understand what you are looking for, so it may not be quite finished yet.

partial-solution.mw

Maple is attempting to produce an answer to an ill-posed problem.  Neither answer satisfies your initial condition D(y)(0) = 0.  You can see that by plotting them:

plot((1/3)*sin(3*t), t=0..2);
plot((1/3*(Heaviside(t)-1/2))*sin(3*t), t=0..2);

That's not a fault of Maple's. Specifying a delta function acting at t=0 is not meaningful. Consider applying it at a later time, as in:

eq := diff(y(t),t,t) + 9*y(t) = Dirac(t-1);

Remove the assignments x:=... and  y:=....  Instead, plot the result through:

odeplot(sol, [u(s)*sin(varphi(s)), -u(s)*cos(varphi(s))], s=0..10);