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

You can set defaut options for all 3-D plots with the command

plots:-setoptions3d('axesfont'= ["TimesNewRoman", 24], 'labelfont'= ["TimesNewRoman", 24]);

These options can be overridden, if desired, on any individual plot command simply by including a different version of the option in the individual command.

You can use plots:-setoptions to do the same thing for 2-D plots.

These can be put into an initialization file which is automatically run whenever a new worksheet is loaded or a restart is executed.

One way:

(Int = rcurry(int, 'AllSolutions'))(1/x, x= a..2);

I have verified conclusively that the bondage number of your example graph is 7 by showing that a certain 7-subset of edges is a bondage set (that part is extremely easy), and using multithreaded processing with the Iterator package to verify exhaustively that none of the 439,381,916 edge subsets with less than 7 edges is a bondage set. With careful multithreading (using 8 threads on my computer), I reduced the average real time to verify whether an edge subset is a bondage set to about 17.3 microseconds (equivalent to about 3.4 million subsets checked per minute).

I'll put the main code and the worksheet separately because the main code is in a Code Edit Region within the worksheet.

MinBondageSet:= module()
option `Author: Carl Love <carl.j.love@gmail.com> 2023-Jun-22`;
local
    #Declaring module locals that are modified in multithreaded code as "thread_local"
    #causes the creation of a separate copy for each thread.
    A::thread_local #modified adjacency array after edge removal
;
export
    #All of this module's exports are essentially "static", but they're not declared
    #"static" because this module is not an "object"
    VL, n, V, A0, #VL = vertex labels, n= # of vertices, V={$n}, A0 = list of vertex neighborhoods
    edges, NE,
    DS, #minimum dominating sets of input graph, as a listlist

    #returns the set of minimum dominating sets of G as a listlist:
    MinDominatingSets:= module()
    local 
        N, s,
        Scans:= subs~(
            _s=~ [0, s, s[-1]],
            proc()
            local N0:= eval(N), Ns, v;
                N:= evaln(N); #garbage collection 
                (
                    for s, Ns in eval(N0) do
                    (
                        for v from _s+1 to n do
                            if (N[s,v]:= Ns union A0[v]) = V then [s,v] fi
                        od
                    )
                    od
                )
            end proc
        ),
        ModuleApply:= proc($)
        local Scan, R;
            N:= table([()= {}]);
            for Scan in Scans do if (R:= Scan()) <> () then return [R] fi od;
            [do R:= Scan() until R <> ()]
        end proc
    ;
    end module,

    GetEdges:= ()-> local k; seq['reduce'= `union`](`{}`~(k, select(`>`, A0[k], k)), k= 1..n-1),

    GetGraphData:= proc(G::Graph) #Assign exports specific to this graph:
        VL:= op(3,G);
        n:= nops(VL);                    
        V:= {$n};                             
        A0:= [seq](op(4,G) union~ `{}`~(V));  
        edges:= GetEdges();
        NE:= nops(edges);
        DS:= MinDominatingSets()
    end proc,

    #This proc essentially makes a local copy of the modified graph,
    #but does so in an efficient way:
    RemoveEdges:= proc(E::{set,list}(set), $)
    option threadsafe;
    local A1:= rtable(A0), e;
         for e in E do A1[e[1]] minus= {e[2]}; A1[e[2]] minus= {e[1]} od;
            A:= seq(A1); #This is why A is thread_local
            return
       end proc,

       #returns true iff no dominating set of the original graph is still dominating after
       #edge removal:
       CheckDS:= proc() option threadsafe; local d; andseq(`union`(A[d]) <> V, d= DS) end proc,
       
    #much slower than CheckDS, but potenially useful for some "greedy" algorithms: selects 
    #all of the original dominating sets that are still dominating after edge removal:
    ReduceDS:= ()-> select(d-> `union`(A[d]) = V, DS),

    ModuleApply:= module()
        local 
         k,        #size of edge subsets being checked
            Min, Max,    #min and max values of k to use  
            C,        #combinations Iterator for k-subsets
          st,        #time at start of an iteration
            found?,     #Flag set in one thread to make the others stop: Has a bondage set been found?

            TimeReport:= proc() #per-iteration userinfo timing report and estimate:
            local T:= time['real']() - st;
                if k = Min or T = 0 then return fi;
            userinfo(
                1, 'MinBondageSet', 'NoName',
                sprintf(
                    "Time to check %d-subsets of edges: %9.1f s; "
                    "expected time to check %d-subsets: %9.1f s.",
                    k, T, k+1, T*(NE-k)/(k+1)
                ) 
            ) 
        end proc
      ;
    export
            CheckSplit:= proc(rnk, num) #main multithreaded procedure:
            option threadsafe;
           local (h,g):= ModuleIterator(Object(C, 'rank'= rnk)), J:= g(), B;
                to num while h() do
                    RemoveEdges((B:= edges[[seq](J+~1)]));
                    if CheckDS() and not found? then found?:= true; return B fi
               until found? #can be set true by any thread
            end proc,

            ModuleApply:= proc(
               G::Graph,    
            { #keyword options:
                #Return after fully processing the input graph but 
                #before examining any edge subset:
                    setup_only::truefalse:= false,

                    #Skip the setup phase because the setup of G is already stored
                    #in the module:
                    no_setup::truefalse:= false,
                    
                    #sizes of edge subsets to check for being bondage sets:
                set_sizes::range({posint, identical()}):= .. 
                },
                $
            )
        uses It= Iterator, TT= Threads:-Task;
        local ne, R;        
            if not no_setup then GetGraphData(G) fi;
            if NE = 0 then return {} fi;
            if setup_only then return fi;

            Min:= lhs(set_sizes); if Min=() then Min:= 1 fi;
            Max:= rhs(set_sizes); if Max=() then Max:= NE fi;
            ne:= binomial(NE, Min);
            found?:= false;
            for k from Min to Max do
                    st:= time['real']();
                    C:= It:-RevolvingDoorCombination(NE, k);
                    R:= TT:-Start('passed', 'Tasks'= [CheckSplit, It:-SplitRanks(ne)[]]);
                    if R <> () then return R fi;
                    ne*= (NE-k)/(k+1); #update binomial(NE, k) -> binomial(NE, k+1) 
                    TimeReport()
            od
            end proc
    ;
    end module,

    #utility that lets the user convert vertex-label indices (default output) back
    #to the original vertex labels,
    LabelVertices:= e-> evalindets(e, integer[1..n], j-> VL[j])
;
end module
:

Download BondageNumber.mw

Suggestions:

The circles for r = 1/3 and r = 2/3 serve no didactic purpose, and thus I think that they're detrimental to the pedagogy.

You need rays for theta = Pi/4, 3*Pi/4, etc.

The vertical lines should not extend beyond the circle. If they go outside the circle, it's easy to mistake them as stray marks. Their color should be more visible.

Your series has no free variables. The concept of radius of convergence doesn't make sense without a free variable such that the convergence of the series depends on that free variable being within a certain "radius" of some point. The concept of power series requires a free variable which is raised to "powers".