Done in Maple 2015.1

restart;

A := <a,b,c>:

a := 1: b := 2: c := 3:

B:=convert(A, list);

B;

                   B := [a, b, c]

                     [1, 2, 3]

The code was edited.

Use  Sum  instead of  sum  to prevent premature calculation.

restart;

Squaring:=proc(Expr)

local f, expr, E, n, r;

f:=op(0,Expr); expr:=op(1,Expr);

if f=Sum then

n:=op([2,1],Expr); r:=op([2,2],Expr);

return Sum(expr^2,n=r) + 2*Sum(Sum(subs(n=_m, expr)*expr, _m=op(1,r)..n-1), n=op(1,r)+1..op(2,r)) else Expr^2 fi;

end proc:

Example:

Squaring(Sum(sin(n*x), n=1..N));

value(%);

Another very simple example:

Squaring(Sum(a[k], k=1..5));

value(%);

Remark.  @Mac Dude  your term  sum(sum(E[m]*E[n],m=1..n-1),n=1..N)  is incorrect because should be  m<>n 

The code was edited.

restart;

seq([b,solve({x+b*y=0, b*x-y=10},{x,y})], b=10..20);

The use of strings can solve your problem. This code is easily rewritten as a procedure and used for similar tasks.

a:=2:  b:=4:  c:=4:  L:=["a","b","c"]:  L1:=map(convert, [a,b,c], string):

S:="a + b*c":  S1:=parse(S): T:=S:

for i to 3 do

T:=StringTools[SubstituteAll](T, L[i], L1[i]);

od:

convert(cat(S, " = ", T, " = ", S1), symbol);

A:=1.330169822*10^11/(6.552475529*10^9*a+7.554142204*10^9*b);

expand(``(numer(A)/10^9)/``(denom(A)/10^9));

For example:

Explore(eval(y*z+x, {x = x0, y = y0, z = z0}), parameters = [x0 = -10 .. 10, y0 = -10 .. 10, z0 = -10 .. 10]);

The list for f(x,y,z)  if  z=-10..10  (for example if x = 1, y = 2):

[seq(eval(y*z+x, {x = 1, y = 2}), z = -10 .. 10)];

@Rouben Rostamian  Your  z = 1/((1+x)*(1 + sqrt(x /a)*arccot(1/sqrt(a/x))))  does not coincide with the initial   z=1/((1+x)*(1+sqrt(x/a)*arccot(sqrt(a/x)))) . 

It seems Maple can not calculate this integral symbolically  even for concrete  a , but it is easily calculated numerically. This calculation is expedient to issue as a procedure with the parameter  a .

Integral:=a->evalf(Int(1/(1+x)/(1+sqrt(x/a)*arccot(sqrt(a/x))), x=0..infinity)):

We can use this procedure to evaluate the integral for the given numbers a and plot the integral vs a.

Examples:

Integral(0.1), Integral(1), Integral(5);

plot(Integral, 0..10, color=red, thickness=2);

If I understand your problem, here is the procedure to solve it. The last parameter of the procedure is N (by default N=100000) to prevent an infinite loop.

FirstReturn := proc(F::procedure, x::realcons, p::realcons, q::realcons, N::posint:=10^5)

local z, p0, q0, t, u, u0, z0;

z:=x; p0:=evalf(p); q0:=evalf(q);

for t from 1 to N do

u:=z; z:=F(z); u0:=evalf(u); z0:=evalf(z);

if (u0<p0 or u0>q0) and (z0>=p0 and z0<=q0) then return [z0, t]  fi;

od;

end proc:

Example of use. Let's try to solve the following problem. Lengthwise of the unit circle from the point (1,0) we shift by 1 (moving counterclockwise). Find the step in which we will be again in range 0 .. 0.0001 (from the point (1,0)).

FirstReturn(x->x+1-floor((x+1)/(2*Pi))*2*Pi, 0, 0, 0.0001);

                                    [0.0000602, 710]

It is known that if the shift (in this example it is 1) incommensurable with the length of the circle (in the example it is 2*Pi), the resulting set of points is dense on the circle.

The procedure  Dist  checks the points  P1  and  P2  whether they belong to the surface of the cube  x=0..1, y=0..1, z=0..1, and then  (following pagan's geometric method), finds the distance between them.

Dist:=proc(P1::list(numeric),P2::list(numeric))

local IsOnCube, d1, d2, d3, d4;

uses geom3d;

IsOnCube:=proc(x,y,z)

if (x>=0 and x<=1) and (y>=0 and y<=1) and (z>=0 and z<=1) and ((x=0 or x=1) or (y=0 or y=1) or (z=0 or z=1)) then return true else false fi;

end proc;

if IsOnCube(op(P1))=false or IsOnCube(op(P2))=false then error "P1 or P2 don't lie on the cube" fi;

point(Q1,P1);  point(Q2,P2);

if ycoord(Q1)<>0 then

  if xcoord(Q1)=0 then rotation(Q1,Q1,Pi/2,line(l,[1/2,1/2,t],t));     rotation(Q2,Q2,Pi/2,line(l,[1/2,1/2,t],t)) else

  if ycoord(Q1)=1 then rotation(Q1,Q1,Pi,line(l,[1/2,1/2,t],t));     rotation(Q2,Q2,Pi,line(l,[1/2,1/2,t],t)) else

  if xcoord(Q1)=1 then rotation(Q1,Q1,-Pi/2,line(l,[1/2,1/2,t],t));     rotation(Q2,Q2,-Pi/2,line(l,[1/2,1/2,t],t)) else

  if zcoord(Q1)=0 then rotation(Q1,Q1,-Pi/2,line(l,[t,1/2,1/2],t));     rotation(Q2,Q2,-Pi/2,line(l,[t,1/2,1/2],t)) else

  if zcoord(Q1)=1 then rotation(Q1,Q1,Pi/2,line(l,[t,1/2,1/2],t));     rotation(Q2,Q2,Pi/2,line(l,[t,1/2,1/2],t)) fi;fi;fi;fi;fi;

fi;

if ycoord(Q2)=0 then return distance(Q1,Q2) else

  if xcoord(Q2)=0 then return distance(Q1,rotation(Q2,Q2,Pi/2,line(l,[0,0,t],t))) else

  if zcoord(Q2)=0 then return distance(Q1,rotation(Q2,Q2,-Pi/2,line(l,[t,0,0],t))) else

  if xcoord(Q2)=1 then return distance(Q1,rotation(Q2,Q2,-Pi/2,line(l,[1,0,t],t))) else

  if zcoord(Q2)=1 then return distance(Q1,rotation(Q2,Q2,Pi/2,line(l,[t,0,1],t))) fi;fi;fi;fi;fi;

if ycoord(Q2)=1 then

    d1:=distance(Q1,rotation(Q2,rotation(Q2,Q2,Pi/2,line(l,[0,1,t],t)),Pi/2,line(l,[0,0,t],t)));

    d2:=distance(Q1,rotation(Q2,rotation(Q2,Q2,-Pi/2,line(l,[t,1,0],t)),-Pi/2,line(l,[t,0,0],t)));

    d3:=distance(Q1,rotation(Q2,rotation(Q2,Q2,-Pi/2,line(l,[1,1,t],t)),-Pi/2,line(l,[0,0,t],t)));

    d4:=distance(Q1,rotation(Q2,rotation(Q2,Q2,Pi/2,line(l,[t,1,1],t)),Pi/2,line(l,[t,0,0],t))); 

return min(d1,d2,d3,d4);

fi;

end proc:

Examples:

Dist([1/2, 1/2, 0], [1, 1, 1]);

                    (1/2)*sqrt(10)

Dist([1/2, 1/2, 0], [1/2, 1/2, 1]);

                                2

Dist([0, 1, 0], [1, 1, 1]);

                           sqrt(2)

Distance_between_points.mw

If you want to have the resulting list has been sorted, then

u:=[x,y,z]; v:=[a,b,x];

[op({op(u), op(v)})];

Matrix(5, RandomTools[Generate](distribution(Normal(10,3)), makeproc=true));

L:=[x1,x2,x3,x4]:

K:=(i::integer) -> `if`(irem(i,4)<>0, L[irem(i,4)], x4);

Examples:

K(0), K(1), K(2), K(3), K(4), K(5);

                  x4, x1, x2, x3, x4, x1

I have simplified and slightly changed your procedure: tr  removed and added the list  L  of colors for faces:

PlotFootball := proc (L::list)

local FootballFaces;

uses geom3d, plottools, plots;

TruncatedIcosahedron(football, point(C, 0, 0, 0), 1);

FootballFaces := seq(polygon(faces(football)[i], color = L[i]), i = 1 .. 32);

display(FootballFaces, axes = none);

end proc:

Example with 4 colors:

PlotFootball([yellow$8, red$8, blue$8, green$8]);

To maintain order, use square brackets instead of curly ones:

plot([sqrt(x+2*sqrt(x-1))+sqrt(x-2*sqrt(x-1)), sqrt(x-2*sqrt(x-1)), sqrt(x+2*sqrt(x-1))], x = -1 .. 5, color=[red,blue,green], thickness=2,numpoints=5000, scaling=constrained);

Green is  sqrt(x-2*sqrt(x-1)) = sqrt(x-1)+1

Blue is  sqrt(x-2*sqrt(x-1)) = abs(sqrt(x-1)-1) = piecewise(x<=2, 1-sqrt(x-1), x>2, sqrt(x-1)-1)

Red = blue + green

I wrote on that Russian forum that the essence of the problem is that the function  g(x)  is the broken line, defined by the points  seq([VE[i], VP[i]], i=1..15)  .  Therefore, the solutions of the equation  g(x)=VP[i]  for each interval  VE[i+1] < x <= VE[i]  coincides with the right boundary of the interval ,  i.e. with  VE[i] . Due to roundoff errors (as Carl wrote) these numbers can be in the following range, so in the range  VE[i+1] < x <= VE[i]  there no solutions.