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;

Obviously, this set is the boundary of the intersection  three regions in the space  -x+y+z <= 1, x-y+z <= 1, and x+y-z <= 1,  that is, the boundary of a polyhedral angle.

A similar example in the plane - the specification of an equilateral triangle by one equation:

plots[implicitplot](max(-y, y-sqrt(3)+sqrt(3)*x, y-sqrt(3)-sqrt(3)*x)=0, x=-1.5..1.5, y=-1..2, thickness=3, gridrefine=3, scaling=constrained);

As for relatively close issues,  see this thread. 

Here is an other solution

ContoursWithLabels(exp(2*x/(x^2+y^2)), -2 .. 2, -2 .. 2, {1, 2, 10, 0.1, 0.5}, [axes = box, color = red, thickness = 2]);

 Edited.  The original function has a singularity at the origin. If you "cut out" it, you can get a pretty good coloring:

f := `if`(x^2+y^2 < 0.015, undefined, exp(2*x/(x^2+y^2))):

ContoursWithLabels(f, -2 .. 2, -2 .. 2, {1, 2, 10, .1, .5}, [axes = box, color = "DarkBlue", thickness = 2], Coloring = [colorstyle = HUE, colorscheme = ["White", "Red"], style = surface]);

restart;

assign(seq(Equ[i]=(u[i](t)=L/2*cos(w*t+phi[i])), i=1..4));

For example:

Equ[1];

                      u[1](t) = 1/2*L*cos(t*w+phi[1])

and so on.

To find all the roots of an analytic function in a certain range it is better to use the  RootFinding[Analytic]  command:

v := unapply(H*sin(w*t), t):

L := 0.080:

H := 0.020:

Vf := 0.3:

w := 10:

zeros1 := sort([RootFinding[Analytic](v(t), t, re = 0.5 .. 2, im = -1 .. 1)]);

zeros2 := sort([RootFinding[Analytic]((D(v))(t), t, re = 0.5 .. 2, im = -1 .. 1)]);

plot([seq([zeros1[i], t, t = -0.01 .. 0.01], i = 1 .. 5), v(t)], t = 0.5 .. 2, color = [yellow$5, blue], linestyle=[3$5,1], thickness=[2$5,1]);

Similarly, you can build  vertical lines at the points of extrema.

in Maple 2015.1

z := x+I*y:

Expr := evalc(abs((1/2)*z^2+3*z-1));

plots[inequal](Expr <= 1, x = 0 .. 0.8, y = -0.4 .. 0.4, scaling = constrained);

Addition: It is better to paint the inside of the region and its border not by the default but using the appropriate options explicitly. Since in this example, we have an open region, its boundary is shown by a dashed line: 

z := x+I*y:
Expr := evalc(abs((1/2)*z^2+3*z-1)):
plots[inequal](Expr <= 1, x = 0 .. .8, y = -.4 .. .4, optionsfeasible = [color = yellow], optionsclosed = [color = blue, linestyle = 6, thickness = 2], scaling = constrained);