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;

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.