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.

For more expressive picture the lightmodel option and some other options are useful:

plots:-display([plot3d([4*sin(t)*cos(s), 2*sin(t)*sin(s), 2*cos(t)], t = 0 .. Pi, s = 0 .. 2*Pi), plot3d([x, x^2, z], x = -4 .. 4, z = -3 .. 4)], scaling = constrained, style = surface, color = "DimGray", axes = normal, lightmodel = light4, orientation = [40, 40], view = [-4.7 .. 4.7, -2.7 .. 3.7, -3 .. 4.7]);

Just use this simple code for your original expression  Expr :

subs({exp(c*t+d*n-d)=exp(c*t+d*n)*exp(-d), exp(2*c*t+2*d*n-d)=exp(2*c*t+2*d*n)*exp(-d)}, Expr);

An example:

restart;

F := (x,y) -> min(x,y)/max(x,y);

x:=2:  y:=1:

D[1](F)(x, F(x,y))+D[2](F)(x, F(x,y))*D[1](F)(x,y);

Addition: alternatively, you can first evaluate this expression for any  x  and  y , and only then to substitute specific arguments, if necessary. The result is the same:

restart;

F := (x,y) -> min(x,y)/max(x,y);

Expr:=D[1](F)(x, F(x,y))+D[2](F)(x, F(x,y))*D[1](F)(x,y);

eval(Expr, {x=2, y=1});

I do not see any bug. Just Maple not fully solve the problem. Here is a workaround:

restart;

Sol:=rsolve({f(n)=0.5*f(n-1)+0.5*f(n+1), f(0)=1, f(6)=0},f(n));

solve(eval(%,n=0)=1, f(5));

f:=unapply(eval(Sol, f(5)=%), n);

Increase the accuracy of the calculations. For example, if  Digits:=50  everything is all right .

Addition: it's interesting that if you use the  solve  command instead of  Eigenvectors  command, you get the best results:

Digits := 50:
p := unapply(CharacteristicPolynomial(M, x), x):
E := [solve(p(x))];
for k to 8 do p(E[k]) end do;