By setting

Digits:= 50;

and adding option complex to your fsolve command (and no other changes), I get this solution in under 2 seconds:

{A = -2.7553365135418814642586082436429575890825402826031,
B = -0.70285804987973303586180028708027467941012949957141}

Note that this solution is real and didn't require initial specification of intervals.

Using the 3rd-party package DirectSearch (available for free download from the Maple Applications Center), and using its command SolveEquations with option AllSolutions, I get 2 solutions: the same one that I show above, and the one reported by Rouben. I'm not confident that this program, or any numeric method, can guarantee the number of solutions to a system of more than 1 non-polynomial equation, even when all univariate restrictions are meromorphic (as is the case here).

@Rouben Rostamian  The procedure can be shortened like this:

tangent_plane:= P-> [P, (sign*primpart@sort@`.`)((<1,3,5>, <x,y,z>)-~'<P>') = 0]:
tangent_plane~(L);

@Laurenso: For your purposes here, the single command sign is equivalent to your signum@lcoeff, but is faster and more robust. It is designed for this type of operation (exact work on polynomials); signum is not.

Here is my version, which uses lists of integers rather than arrays. It's a very short code. The run time is just slightly longer than dharr's.

restart
:
nextone:= L-> local i:= 2, j, n:= nops(L);
    ["", (do for i from (j:= i) + 1 to n while L[i]=L[j] do od; i-j, L[j] until i > n)]
:
LookAndSay:= (n::posint, start::nonnegint)-> local L:= ["", start];
    <start, (to n-1 do parse(cat((L:= nextone(L))[])) od)>
:
CodeTools:-Usage(LookAndSay(50,1));
memory used=309.93MiB, alloc change=27.45MiB,
cpu time=3.52s, real time=3.52s, gc time=281.25ms

Regarding the timing of your Maple code (the StringTools-based code in your Question): You can't blame the compiled StringTools code. The vast majority of the time is spent doing garbage collection at the Maple kernel level, not in the compiled code. This is due to the large number of temporary substrings that you create.

You just needed some minor syntax corrections:

R:= plots:-implicitplot(
    [x = -3/4, x = -1/3, x = 1/3, x = 3/4, y = -3/4, y = -1/3, y = 1/3, y = 3/4], 
    x = -3 .. 3, y = -3 .. 3, color = red, scaling = constrained, gridrefine = 3
);
z:= x + y*I:
plottools:-transform(unapply(evalc([(Re, Im)(1/z)]), x, y))(R);

The circles don't show as complete circles because the lines are not complete (infinitely long) lines; they're line segments.

To get a single solution, you just need to make 3 trivial changes:

1. Correct beta_1 = 0.2 to beta_1:= 0.2;
2. Correct all pi to Pi;
3. Change solve to fsolve.

If you need all solutions, or all real solutions, it's a little trickier.

The example code that you show

f[N1_]=sum[1/(n^3*sin^2[n]), {n,1,N1}];
DiscretePlot[f[x], {x,0,400},PlotRange->All];

is Mathematica, not Maple.

For Maple, you can do a continuous plot of it as

plot(eval(sum(exp(-lambda)*lambda^n/n!, n= 0..N1), lambda= 15.4), N1= 0..50);

If you insist on a discrete plot, simply replace N1 with trunc(N1) inside the sum:

plot(eval(sum(exp(-lambda)*lambda^n/n!, n= 0..trunc(N1)), lambda= 15.4), N1= 0..50);

You need to set your working directory to something that your OS gives you permission to write a file to. This can be done with the currentdir command. For example, on my computer, this works:

currentdir("/users/carlj/desktop");
ExportMatrix(matlabData, A, target= MATLAB, format= rectangular, mode= ascii);

Let's suppose that:

1. A is any of the arguments of exp that occur in your expressions.
2. You're able to write a procedure TransA that takes any such A and returns what it should be replaced with.
3. Expr is any expression or list, set, etc., of expressions containing some exp functions that you want to change.

Then, this command will do what you want:

subsindets(Expr, specfunc(exp), exp@TransA@op)

Your procedure TransA must return something for any A that it's passed. If it simply returns an A unchanged, that's fine.

If A itself contains exp subexpressions, there's no need for TransA to transform their arguments, because the subsindets will always process them in order from most deeply nested to least deeply nested.

There's nothing special about exp here. The above process will work for any function of a single argument. With slight modifications, it'll also work for multi-argument functions.

Since the probability is very high that the desired condition will be met anyway, I think that the best strategy is to check that the condition is met, and, if it's not, regenerate the matrix. Like this (using 1-D input):

do
    A:= LinearAlgebra:-RandomMatrix((10,10), 'generator'= -10..10)
until andseq(rtable_scanblock(A, [i, ..], 'NonZeros') > 2, i= 1..10);

You can also generate the rows satisfying the condition one at a time and catenate them, like this (1-D input) (I've added the density= 0.3 to give the problem a more-realistic level of difficulty):

A:= <(
    to 10 do
        do
            R:= LinearAlgebra:-RandomVector(
                10, 'generator'= -10..10, 'density'= 0.3
            )^+
        until rtable_scanblock(R, [], 'NonZeros') > 2
    od
)>;

You wrote:

• The bases that exit for one selected period (your proc) are pairwise conected to another base in the set.  I can see this in the plots.

  So  it's possible to run only the half of the loop with some extra calculations for pairfinding. I think this may reduce the time to 60% of running the fulll loop.

Yes, both of those observations are essentially correct, but actually the time reduction is usually much greater. The "orbits" of a given base with respect to a given period are not just pairs but have size totient(period). So, for example, for period = 11, the bases come in groups of size totient(11) = 10. Knowing any one base in that orbit, the other 9 are trivial to generate. Here is a program for it:

restart:
Periods:= module()
export
    #prime factorization in convenient list form (see ?ifactors):
    Ifactors:= (n::posint)-> (thisproc(n):= ifactors(n)[2]),

    #prefactored form of the totient, making it more efficient to compute the 
    #Carmichael number at the same time:
    TotSeq:= (n::posint)-> 
    local p; 
        (thisproc(n):= [seq]((p[1]-1)*p[1]^(p[2]-1), p= Ifactors(n))),

    #Euler's totient function:
    Totient:= (n::posint)-> (thisproc(n):= mul(TotSeq(n))),

    #Carmichael's lambda function:
    Carmichael:= (n::And(posint, Not(1)))-> 
    local T:= TotSeq(n);
        (thisproc(n):= ilcm(`if`(n::even and T[1]>2, [T[1]/2, T[2..][]], T)[])),

    #list of prime divisors of n:
    PrimeDivisors:= (n::posint)-> (thisproc(n):= op~(1, Ifactors(n))),

    #Carmichael number reduced by its prime divisors:
    PreCarmichael:= (n::And(posint, Not(1)))->
    local L:= Carmichael(n);
        (thisproc(n):= iquo~(L, PrimeDivisors(L))),

    (* Find any element of order Carmichael(n). The optional table Ex contains
    elements to NOT check. The optional 3rd argument St is the element 
    to start with. *) 
    Generator:= proc(
        n::And(posint, Not({1,2})), Ex::table:= table('sparse'), St::posint:= 2
    )
    local E:= PreCarmichael(n), a, e;
        for a from St to n-1 do
            if Ex[a]=0 and igcd(a,n)=1 and andseq(a&^e mod n <> 1, e= E) then 
                return a 
            fi
        od;
        FAIL
    end proc,

    #Find the elements of order ord. The optional 3rd argument NG is the
    #maximum number of generators to use.
    Bases:= proc(
        n::And(posint, Not({1,2})),
        ord::And(posint, Not(1)),
        NG::posint:= infinity
    )
    local 
        Ex:= table('sparse'), L:= Carmichael(n), LP:= iquo(L, ord), 
        E:= select(igcd=1, [$1..ord-1], ord), EG:= select(igcd=1, [$1..L-1], L), 
        G, e, g:= 1, r, ng:= 0, h, TT:= Totient(Totient(n))
    ;
        if irem(L, ord) <> 0 then return {} fi;
        E:= E[2..] -~ E[..-2];
        EG:= EG[2..] -~ EG[..-2];
        {
            to NG while ng <> TT and (g:= Generator(n, Ex, g+1)) <> FAIL do
                G:= e-> (thisproc(e):= g&^e mod n);
                Ex[(h:= g)]:= 1; ng++;           
                for e in EG do 
                    if Ex[(h:= h*G(e) mod n)]=0 then Ex[h]:= 1; ng++ fi 
                od;
                r:= g&^LP mod n;
                G:= e-> (thisproc(e):= r&^e mod n); 
                ((h:= r), seq((h:= h*G(e) mod n), e= E))
            od
        }
    end proc  
;
end module
:
den:= 11842585:  #same number that you used
NumberTheory:-Divisors(Periods:-Carmichael(den));
{1, 2, 3, 4, 6, 7, 11, 12, 13, 14, 21, 22, 26, 28, 33, 39, 42, 
  44, 52, 66, 77, 78, 84, 91, 132, 143, 154, 156, 182, 231, 273, 
  286, 308, 364, 429, 462, 546, 572, 858, 924, 1001, 1092, 1716, 
  2002, 3003, 4004, 6006, 12012}

bases11:= CodeTools:-Usage(Periods:-Bases(den, 11));
memory used=287.81MiB, alloc change=0 bytes, 
cpu time=4.28s, real time=3.79s, gc time=703.12ms

bases11 := {207496, 299716, 338141, 484156, 514896, 561006, 
  568691, 737761, 853036, 1029791, 1098956, 1237286, 1367931, 
  1452466, 1544686, 1590796, 1598481, 1752181, 1767551, 1844401, 
  2105691, 2128746, 2159486, 2267076, 2282446, 2359296, 2482256, 
  2620586, 2781971, 2797341, 3135481, 3143166, 3189276, 3312236, 
  3389086, 3573526, 3658061, 3673431, 4026941, 4188326, 4418876, 
  4457301, 4541836, 4603316, 4634056, 4733961, 4933771, 5202746, 
  5248856, 5448666, 5487091, 5663846, 5717641, 5733011, 5871341, 
  5963561, 6086521, 6148001, 6224851, 6247906, 6401606, 6601416, 
  6662896, 6747431, 6793541, 6901131, 7031776, 7208531, 7254641, 
  7277696, 7308436, 7431396, 7546671, 7631206, 7723426, 7777221, 
  8207581, 8292116, 8307486, 8338226, 8445816, 8660996, 8722476, 
  8799326, 8807011, 8822381, 8960711, 8976081, 9052931, 9091356, 
  9175891, 9268111, 9314221, 9321906, 9337276, 9368016, 9490976, 
  9606251, 9690786, 9752266, 9783006, 9852171, 9882911, 10082721, 
  10205681, 10397806, 10505396, 10636041, 10782056, 10812796, 
  10858906, 10912701, 11020291, 11035661, 11296951, 11373801, 
  11550556, 11627406, 11665831, 11811846}

The last-used index of a for loop is still accessible after the loop ends. At that point, you can, for example, subtract 1 from it and use it. Like this:

for i to nops(testseq) while testseq[i] < 15 do od:
const:= `if`(i=1, FAIL, testseq[i-1]);

If all elements satisfy the while condition, then i will equal nops(testeq) + 1 at the end, and the last command will still work.

Yes, it can be done as in this example:

for j to 99 do
    p:= ithprime(j);
    try
        timelimit(1, ifactor(2^p-1))
    catch "time expired":
        printf("Time expired at j=%d; proceeding to next j.\n", j);
        next  #use next j
    end try
od:

However, note that there are some commands that may use far more than their allowed amount of time before the timelimit command "catches" and stops them.

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.