It's not mentioned on its help page, but sequences of plots structures are also accepted by plots:-display. These may be followed by additional options.

Many Maple commands accept additional forms of input that they aren't documented to accept. I'm not quite sure if that's a good thing or a bad thing. 

The code for the module named action is in a part of the worksheet called "Startup Code", as mentioned by Joe Riel. When you load the worksheet, you need to grant permission to Maple to execute that code. You'll get a popup dialog box:

• "Autoexecute Warning: This worksheet contains content that will execute automatically. Do you wish to proceed?"

You need to click "Yes".

There are always a few language additions in each new version of Maple. In the vast majority of cases, older code continues to work. (That is a design goal called "backward compatibility".) Once I grant permission, it seems to work. Let me know if you have any more problems.

Although you haven't said this explicitly, the way that you present the Question suggests that you expect to count those black points via some sort of image processing of the final plot. However, it's much easier to just save those points when they are found by your numeric algorithm. They are, of course, those points for which the algorithm returns 0.

There are two ways to do this. The one you use will be determined by whatever after-processing you want to do. They are both very easy to implement, but the 1st is a little bit easier than the 2nd: 

1. Simply count those points.
2. Save those points.

Both cases involve adding exactly one line of code to Hafiz1Count, immediately before the return 0, and exactly one line outside the procedure before the plot3d is called. To simply count them, the line in the procedure is
    :-BlackCt++;   #or :-BlackCt:= :-BlackCt+1
and the outside line is
    BlackCt:= 0;

To save the points, the line inside the procedure is
    :-BlackZ,= x0+y0*I;  # Yes, the operator is "comma equal"!
the outside line is
    BlackZ:= Array(1..0, datatype= complex[8]);
and the number that have been saved can be retrieved by
    numelems(BlackZ);

Both of my suggestions for the inside line require that the procedure be entered in 1D input.

In both cases, the outside line, which initializes the counter or array, must be re-executed before redoing any plots that use Hafiz1Count. 

You wrote:

•  I would expect the way to do this would be ans_all[1].

That would be true except for cases where allvalues only returns 1 solution, which is what is happening in this case. If you want an expression that'll work regardless of whether there's 1 or more than 1 solution, use

[ans_all][1]; 

Your inner procedure computus does not use any lexical variables. Thus, it is utterly trivial to de-nest it! With a few trivial syntactic (not algorithmic) modifications, it can be compiled. 

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]
);