The basic technique for creating dynamic nested loops is to apply the foldl (fold-left) command to seq or one of its companion looping commands such as add. This is what @sursumCorda did. As he showed, making the loops dynamic by folding doesn't make them any faster, nor should you expect it to. But I did find 2 major ways to make the code faster, speeding up yours and @sursumCorda 's by a factor of at least 40 (i.e., 40 times faster).

The first improvement was to eliminate your innermost loop. That loop is simply counting (by repeatedly adding 1) a sequence of consecutive integers. If you know the integer bounds, the count can be done by simple subtraction of those bounds, plus 1. That's the significance of this expression, which is the innermost addend in my add loops:

max(0, min(_v||n1 + 1, V[-1]) - Ceil(V[-1]/2))

The second improvement is due to a Maple idiosyncracy: The command trunc is many times faster than ceil (or floor or the other commands of that ilk). By making a slight adjustment to trunc, it can replace ceil:

ceil(x) = trunc(x) + `if`(trunc(x)=x or x < 0, 0, 1)

With these changes, my final code is able to do the n=6 case in under 9 seconds.

nequaln:= (n::And(posint, Not(1)))->
local 
    i, %add, n1:= n-1, 
    V:= (n-2)*24 -~ (0, seq['scan'= `+`](_v||i, i= 2..n1)), #seq of partial sums

    #Maple's library 'ceil' has a highly symbolic aspect that slows it significantly. (You
    #can confirm this by showstat(ceil).) But 'trunc' is kernel builtin and doesn't have
    #this problem. Thus, I wrote this version of 'ceil' that only uses 'trunc'.
    Ceil:= x-> local T:= trunc(x);  
        `if`(T::integer, (thisproc(x):= T + `if`(T=x or x<0, 0, 1)), 'procname'(x))
;
    (eval@subs)(
        _v||1= V[1] - n1, %add= add,
        foldl(
            %add,
            max(0, min(_v||n1 + 1, V[-1]) - Ceil(V[-1]/2)),
            seq(_v||i= Ceil(V[i-1]/(n-i+2)).._v||(i-1), i= n1..2, -1)
        )        
    )           
:
[seq](CodeTools:-Usage(nequaln(k)), k= 2..6);
memory used=14.30KiB, alloc change=0 bytes,       #Usage for k=2
cpu time=0ns, real time=2.00ms, gc time=0ns

memory used=285.20KiB, alloc change=0 bytes,      #Usage for k=3
cpu time=15.00ms, real time=14.00ms, gc time=0ns

memory used=475.34KiB, alloc change=0 bytes,      #Usage for k=4
cpu time=0ns, real time=9.00ms, gc time=0ns

memory used=13.45MiB, alloc change=0 bytes,       #Usage for k=5
cpu time=79.00ms, real time=80.00ms, gc time=0ns

memory used=1.58GiB, alloc change=0 bytes,        #Usage for k=6
cpu time=8.58s, real time=8.35s, gc time=437.50ms

                  [0, 48, 816, 10642, 117788]

You're trying to use a variable named delta[h]_range. That's not allowed as an unquoted  variable name. But, if you wish, you can use back-quotes (aka left-quotes) to make anything at all a variable name:

`delta[h]_range`:= 0.55 .. 0.75;

Click anywhere in the output field--no need to highlight it--and press Ctrl-Delete.

It can be done like this:

period:= 6:  x:= 78:
d:= period /~ NumberTheory:-PrimeFactors(period):
select(
    b-> igcd(x,b)=1 and b&^period mod x = 1 and andseq(b&^i mod x <> 1, i= d),
    [$2..x-1]        
);
                    [17, 23, 29, 35, 43, 49]

The asymptotic time complexity with respect to period can be further reduced by using some kind of "factor tree" to compute the powers b &^ i mod x, where i runs over the proper divisors of period. There's no benefit to this if the period is square-free, as is 6. 

Edit: The idea mentioned in the previous paragraph is now implemented by the new 2nd line above:

d:= period /~ NumberTheory:-PrimeFactors(period);

You can put the evaluation of the Newton iteration in a try-catch statement to trap both singularities of f and points where its derivative is 0. 

NM2:= proc(
    f, xold::complexcons, {precision::positive:= 10.^(1-Digits), max_iters::posint:= 9}
)
local x0:= xold, x1:= x0+2*precision, Df:= D(f), k;
    for k to max_iters do
        try 
            (x0,x1):= (x1, evalf(x1 - f(x1)/Df(x1)))
        catch:
            error "Singularity at %1", x1
        end try;
        userinfo(1, procname, 'NoName', x1)
    until abs(x0-x1) < precision; 
    if k > max_iters then error "Maximum iterations exceeded" fi;
    x1 
end proc
:
p:= randpoly(x, degree= 9, dense);
            9       8       7       6       5       4       3
   p := 72 x  + 37 x  - 23 x  + 87 x  + 44 x  + 29 x  + 98 x 

            2            
      - 23 x  + 10 x - 61

f:= unapply(p, x):
infolevel[NM2]:= 1:
r:= NM2(f, 1.3);
1.121934407
.9609516325
.8274608051
.7423627912
.7136509949
.7110238945
.7110039767
.7110039755
.7110039755
                       r := 0.7110039755
f(r);
                               0.
r:= NM2(f, 1.+I);
.9145078329+.9127717630*I
.8599212591+.8460365442*I
.8403634004+.8059364156*I
.8411389647+.7949972141*I
.8416536274+.7947359908*I
.8416536540+.7947375421*I
.8416536539+.7947375422*I
               r := 0.8416536539 + 0.7947375422 I
f(r);
                                    -7  
                        0. + 1.42 10   I

I think that the following is what you have in mind, although it's not quite correct to say that it's the surface common to the two cylinders. Rather, it's the portion of each cylinder that's inside the other.

plots:-display(
    plot3d(  # y^2 + z^2 = 1:
        [[r, t, sqrt(1-r^2*sin(t)^2)], [r, t, -sqrt(1-r^2*sin(t)^2)]],
        r= 0..1, t= -Pi..Pi, coords= cylindrical
    ),
    plot3d(  # x^2 + y^2 = 1:
        [1, t, z], 
        t= -Pi..Pi, z= -sqrt(1-sin(t)^2)..sqrt(1-sin(t)^2), coords= cylindrical
    ),  
    scaling= constrained, labels= [x,y,z]
);

You seem to be having trouble distinguishing a procedure's input from its output. You want a procedure whose input is the integer dimensions of the matrix and whose output is the matrix itself. 

Matrix01:= proc(m::nonnegint, n::nonnegint:= m)
    Matrix((m,n), (i,j)-> modp(i+j, 2)) 
end proc
:
#Examples:
Matrix01(4);
Matrix01(3,5);

A fundamental motivation of the discipline of Combinatorics is that counting combinatorial structures can often be done using far less computational effort than generating the structures. Counting partitions is like that.

It can be done quite efficiently with a very short recursive procedure:

restart:
CountSums:= (k::posint, N::set(posint))->
    ((k,n)-> (thisproc(args):= 
        `if`(k*n > 0, thisproc(k, n-1) + thisproc(k-N[n], n), `if`(k=0, 1, 0))
    ))(k, nops(N))
:  
CodeTools:-Usage(CountSums(1000, {2,3,6,9}));
memory used=421.16KiB, alloc change=0 bytes,
cpu time=15.00ms, real time=7.00ms, gc time=0ns

                             526848
CountSums(9, {2,3,6,9});
                               4
CountSums(8, {2,3,6,9});
                               3
#Verify correctness:
CodeTools:-Usage(combinat:-numbpart(100));
memory used=2.14MiB, alloc change=0 bytes,
cpu time=31.00ms, real time=33.00ms, gc time=0ns

                           190569292

CodeTools:-Usage(CountSums(100, {$100}));
memory used=1.18MiB, alloc change=0 bytes, 
cpu time=16.00ms, real time=13.00ms, gc time=0ns

                           190569292

So, for this example, my procedure's time beats the library procedure combinat:-numbpart.

@Pepini You wrote:

• I thought about use parametrization of torus:
  plot3d([(R + r*cos(t)*cos(s), (R + r*cos(t)*sin(s), r*sin(t)], t = 0 .. 2*Pi, s = 0 .. 2*Pi)

You are missing some parentheses in that parametrization: The R also needs to be multiplied by the cos(s) or sin(s). That omission may be why your plot wasn't satisfactory as a surface of revolution. Here's a correction:

r:=  t-> 1 + sin(8 + sin(16 + sin(32*t)))/4:  #This works for any polar curve in form r = f(t)
R:= 2:
plot3d(
     [(R + (r*cos)(t))*~(cos,sin)(s), (r*sin)(t)], t= -Pi..Pi, s= -Pi..Pi,
     scaling= constrained,  #usually a better option than using 'view'
     grid= [101,49]   #Increase points used for t to show wrinkling
);

This is a figure homeomorphic to an ordinary torus, not the "flat torus" that @sand15 gave links to a PDF about and you showed a picture of. And it has only one level and direction of what those authors call "corrugation" (you might call it "wrinkling"). I think that 2 directions would be the minimum for a reasonable plot, but 1 direction is the most that can be accommodated by a surface of revolution. 

Here's a Matrix with all the information you need. You could display it more elegantly with DocumentTools, I suppose.

Download NormalTable.mw

Perhaps you will understand better if things are made more explicit rather than using the geometry package. MaplePrimes won't let me display this worksheet inline at the moment, but you can download it...

Download Ellipse.mw

...or just copy and paste this code into a Maple session:

Dist:= (P,Q)-> sqrt(add((P-Q)^~2)):  #Euclidean distance between points P and Q
Mid:= (P,Q)->  (P+Q)/2:              #Midpoint between P and Q

P:= <x,y>:                       #An arbitrary symbolic point on the ellipse
F:= [<0,-5>, <5, 0>]:  MA:= 12:  #Foci and major-axis length

(* An ellipse's defining property: The sum of the distances from any of its
points to its foci is constant. I'll call that constant C for the moment. 
We'll use that other piece of given information, the length of the major axis,
to figure out C in a moment. *)

Ell_eqn:= C = add(map(Dist, F, P));

# Perhaps you'll accept without further explanation that the line through 
# the foci (which is collinear with the major axis) is
MA_line:= y = x - 5;

(* The next step to getting C is finding numerically any point on the ellipse. 
Two such points are the points on that line that are at distance MA/2 from the
midpoint of the foci. We solve for them: *)

%solve({MA_line, Dist(P, Mid(F[])) = MA/2}, explicit);
sols:= value(%); 

# As suspected, there are 2 solutions. We only need 1, and I arbitrarily choose
# the first. To find C, we just use those x and y in the ellipse's equation:

eval(Ell_eqn, sols[1]);
eval(Ell_eqn, %):  #Put that in for C

# Subtract one side of the equation from the other to make an expression 
# implicitly equated to 0:

(lhs-rhs)(%);    
evala(Norm(%));      #Convert that to a polynomial
(primpart/sign)(%);  #Cancel common integer factors of the coefficients
#Final:
             119*x^2 - 50*x*y + 119*y^2 - 720*x + 720*y - 1584 

Your expression involves extreme underflows and overflows that can't be handled by the default hardware floats. So, the easy correction is

UseHardwareFloats:= false: 

An alternative is to simplify the expression to

f2:= exp(5002.646494 + 2499*ln(t) - 5000.*sqrt(t));

Then hardware floats can be used by plot.

My preferred syntax for situations like this is decisions[`if`(1 > T, 1, 2)].

I told you before that you need to change Gamma to GAMMA. You said that you did, but I still see Gamma.

Try entering it as A[b,c]. This makes it an indexed variable. The output will be subscripted like you want. An indexed variable is not quite the same thing as the atomic variable that you were trying to create, but it'll work if you do not assign values to A, b, or c.

The characters that you see are not apostrophes. Note that they lean to the left. They're called "back quotes" or "accents graves" or, in Maple, "name quotes". They let you make any string of characters, no matter how weird, into a variable.

However, you shouldn't be seeing those back quotes in your output. If I enter 

`A__b,c`;

then my output appears without quotes as A_b,c