1) See  LinearAlgebra[HilbertMatrix]  command.

2) Your matrix  (2*2 matrix which its transpose is egual to its inverse)  is a rotation matrix or composition of a rotation and a reflection:  

Matrix([[cos(alpha), -sin(alpha)], [sin(alpha), cos(alpha)]]);  # A rotation matrix

Matrix([[1,0], [0,-1]]) . Matrix([[cos(alpha), -sin(alpha)], [sin(alpha), cos(alpha)]]);  # Composition of a rotation and a reflection

We have a system with 27 unknowns and 54 variables, that is 27 variables are parameters. If you give some values to these parameters, fsolve  easy solves the system. Replace the last line of your code by

ABC := fsolve(eval(map(t->Equate(op(t))[], [eq2 = AA1, eq3 = AA2, eq4 = B2]), {aaa1 = 1, aaa2 = 2, aaa3 = 3, aaa4 = 4, aaa5 = 5, aaa6 = 6, aaa7 = 7, aaa8 = 8, aaa9 = 9, aab1 = 10, aab2 = 11, aab3 = 12, aab4 = 13, aab5 = 14, aab6 = 15, aab7 = 16, aab8 = 17, aab9 = 18, aabb1 = 19, aabb2 = 20, aabb3 = 21, aabb4 = 22, aabb5 = 23, aabb6 = 24, aabb7 = 25, aabb8 = 26, aabb9 = 27}));

{a1 = 1.461172682, a2 = 0.4611726819, a3 =0 .4611726819, a4 = 8.077654636, a5 = 9.077654636, a6 = 8.077654636, a7 = -8.538827318, a8 = -8.538827318, a9 = -7.538827318, b1 = -3.221842064, b2 = -2.273083473, b3 = -2.324324882, b4 = 10.46946092, b5 = 8.571943737, b6 = 8.674426555, b7 = -7.673550001, b8 = -6.724791410, b9 = -6.776032819, c1 = 31.22184206, c2 = 31.27308347, c3 = 32.32432488, c4 = 20.53053908, c5 = 23.42805626, c6 = 24.32557344, c7 = 41.67355000, c8 = 41.72479141, c9 = 42.77603282}

This is an answer to the second question. The answer to the first question is obvious.

Basic variables corresponds to the maximum set of linearly independent columns of the leading matrix of the system. The remaining variables are the free variables, which we consider as parameters. Selection of basic variables is ambiguous. I think the default Maple selects the first set in lexicographical order. For example in solving of  {x-y + z = 2, 2 * x + y-z = 3}  Maple expresses  x  and  y  through  z . But for the system  {x-y+z=2, 2*x-2*y-3*z=5}  x  and  z  are expressed through  y , because the first 2 columns are linearly dependent. 

You can choose yourself the basic variables (if the corresponding columns are linearly independent) as in the example:

solve({x-y+z=2, 2*x+y-3*z=5}, {x,z});

                                   {x = 11/5+2/5*y, z = -1/5+3/5*y} 

Maple does not directly plot the set. Some simplification is necessary to this end:

evalc(eval(abs(z+1+I)>=abs(z-1-I), z=x+I*y));

map(t->t^2, %);

plots[inequal](%, x=-5..5, y=-5..5, color="LightGreen");

twinprimes:=proc(a,Lim)

local a0, b;

a0:=`if`(a::prime, a, nextprime(a));  b:=nextprime(a0);

while b<=Lim do

if b-a0=2 then  lprint([a0,b])  fi;

a0:=b; b:=nextprime(a0);

od;

NULL;

end proc:

Example of use:

twinprimes(10,1000);

Edited.

Consider the graphical solution, using the continuity of the function  (a,x)->log[2*abs(x-a)](abs(x+a)+abs(x-a))-1  in some regions. To find these regions we consider the lines on which no solutions:

log[2*abs(x-a)](abs(x+a)+abs(x-a)) = 1  is equivalent to  abs(x+a) = abs(x-a)

2*abs(x-a) = 0  is equivalent to  x=a

2*abs(x-a) = 1  is equivalent to  x=a+1/2  or  x=a-1/2

abs(x+a)+abs(x-a)=0  is equivalent to  one point  a=0  and  x=0

The plot of these lines:

 plots[implicitplot]([abs(x+a)=abs(x-a), x-a=0, 2*abs(x-a)=1], a=-1..1, x=-1..1, color=[red,blue,green], thickness=4, gridrefine=5, axes=normal);

It is enough to consider only the regions for  a>0  because the mapping  (a,x)->(-a,-x)  does not change the inequality. 

Consider one point in each of the 6 regions and the truthfulness of the inequality in these regions:

R1, R2, R3, R4, R5, R6 := [0.5,-0.5], [0.2,-0.2], [0.5,0.2], [1,0.2], [0.5,0.7], [0.5,2]:

seq(is(eval(log[2*abs(x-a)](abs(x+a)+abs(x-a))<1, [a,x]=~(R||i))), i=1..6) ;

                                  true, false, true, false, true, false

So the final solution:

If  a<=-1/2  then  x>0  or  a<x<a+1/2  or  a-1/2<x<a

If  -1/2<a<0  then  x>a+1/2  or  a<x<0  or  a-1/2<x<a

If  0<a<=1/2  then  x<a-1/2  or  0<x<a  or  a<x<a+1/2

If  a>1/2  then  x<0  or  a-1/2<x<a  or  a<x<a+1/2

Addition: visualization of all the solutions (yellow regions):

P:=plots[implicitplot]([abs(x+a)=abs(x-a), x-a=0, 2*abs(x-a)=1], a=-1..1, x=-1..1, color=[red,blue,green], thickness=2, gridrefine=3, axes=normal):

A:=plots[inequal](x>piecewise(a<-1/2,0,a>-1/2 and a<0,a+1/2), a=-1..0, x=0..1, color=yellow):

B:=plots[inequal]({x<piecewise(a<-1/2,a+1/2,a>-1/2 and a<0,0), x>a-1/2}, a=-1..0, x=-1..0, color=yellow):

C:=plots[inequal]({x>piecewise(a>0 and a<1/2,0,a>1/2,a-1/2),x<a+1/2}, a=0..1, x=0..1, color=yellow):

E:=plots[inequal](x<piecewise(a>0 and a<1/2,a-1/2,a>1/2,0), a=0..1, x=-1..0, color=yellow):

plots[display](P,A,B,C,E);

Clarification: on the blue line there are no solutions. 

The following code finds all the fours of integer points on your sphere (for which all 12 coordinates are different) and which form a square:

restart;

n := 0:

for x from -13 to 17 do

for y from -11 to 19 do

for z from -9 to 21 do

if (x-2)^2+(y-4)^2+(z-6)^2 = 225 then n := n+1; S[n] := [x, y, z] fi;

od: od: od:

S := convert(S, set):

S1 := combinat[choose](S, 3):

n:=0:

for s in S1 do

P:=combinat[permute](s);

for p in P do

if add((p[2]-p[1])*~(p[3]-p[1]))=0 and add((p[2]-p[1])^~2)=add((p[3]-p[1])^~2) then n:=n+1;

L[n]:={op(p),p[2]+p[3]-p[1]}; break fi;

od: od:

L:=convert(L, set):

Squares:=select(p->nops(op~(p))=12, L); nops(%);

 Squares_on_sphere.mw

restart;

k:=x->piecewise(type(x,even),x/2, 3*x+1);

for m from 1 to 60 do

s:='s';  m1:=m;  n:=1;  s[1]:=m1;

while m1>1 do n:=n+1;  m1:=k(m1); s[n]:=m1;  od:

print(op(convert(s, list)));

od: 

Addition: If you need to check the conjecture for any set of numbers, these numbers are convenient to specify as a list or a set.

Example for 2 numbers  10^6  and  10^9 :

restart;

k:=x->`if`(type(x,even),x/2, 3*x+1):

for m in [10^6, 10^9] do

s:='s';  m1:=m;  n:=1;  s[1]:=m1;

while m1>1 do n:=n+1;  m1:=k(m1); s[n]:=m1; od:

print(op(convert(s,list)));

od:

1) If we use a sequence (or a list of some more structures) in the end of the procedure, it is possible to do without return.

2) A colon is equivalent to a semicolon inside the body of a procedure.

3) Semicolon or colon may not be used before  end  (or  end proc).

LTTS:=proc(ff)

           local ll,r,r1,r2,r3;

           ll:=rhs(ff)-lhs(ff);

                     solve01(ll),

                     solve02(ll),

                     solve03(ll),

                     solve04(ll),

                     solve05(ll)

          end:

See help on  PolyhedralSets[Volume]  and  PolyhedralSets[Area]  commands. It is in Maple 2015 only.

For other versions Maple there are commands in the plane only:  simplex[convexhull]  and  geometry[convexhull]

Addition: an example

with(PolyhedralSets):

Solid := PolyhedralSet([[0, 0, 0], [0, 1, 0], [1, 1, 0], [0, 1, 1]]);

Volume(Solid);

Plot(Solid, axes = normal);

Without any assumptions only by scaling:

E:= 3^(-(1/2)*n)*2^((1/6)*n)-2^((2/3)*n)*6^(-(1/2)*n):

simplify(subs(n=6*k, E));

                                                 0

Addition:  Since the mapping  z -> 6*z  is a bijection  C  on  C  then the identity  E=0  is true for any complex number  n .

zahl:=1234567:

(floor(sqrt(zahl))+1)^2;

                                         1236544

With a loop:

zahl:=1234567:

 for n from zahl+1 while not type(sqrt(n), integer)  do 

end do:

n;

                                         1236544

Edited.

@litun 

See the corrected variant  Roots.mw  of your file  test.mw

All the roots (real and imaginary) of an analytic function  in specific ranges you can get by  RootFinding[Analytic]  command.

An example (parameters  B .. beta  are taken arbitrarily):

f:=(A,B,Psi,K1,K2,K3,alpha,beta,m) -> A-B*((m+1)/Psi+(K1*(1-beta^(m+1))+K2*alpha^(m-1))/K3/alpha^(m-1));

for A from 0 to 1 by 0.1 do

RootFinding[Analytic](f(A,1,2,3,4,5,1,2,m), m, re=-10..10, im=-10..10);

od;

We see that for each of the first 9 values of parameter  A  the function has 2 real roots, for A = 0.9  and  A=1  - imaginary roots. 

assuming  option will help you:

int(sigma*exp(-sigma*x)/(1+exp(-sigma*x))^2, x = -infinity .. infinity)  assuming sigma > 0;

int(exp(-x/sigma)/(sigma*(1+exp(-x/sigma))^2), x = -infinity .. infinity)  assuming sigma > 0;

int(exp(-(x-mu)/sigma)/(sigma*(1+exp(-(x-mu)/sigma))^2), x = -infinity .. infinity)  assuming sigma > 0,  mu::realcons ;

                                                                            1

                                                                            1

                                                                            1

Edited.