If you read through the code showstat(ifactor), you'll see that the specification of a method in the second argument is open ended. In other words, any name meth can be given, any then ifactor will pass your number on to `ifactor/meth`. So, you're not restricted to the methods mentioned at ?ifactor. As I said earlier, the Maple 2021 library has 30 procedures so named, and the tags meth can be extracted like this:

Meths:= cat~(``, index~(LibraryTools:-PrefixMatch(`ifactor/`), 9..-3));
[PollardP1, PollardP1/ext, PollardP1/tree, PollardRhoBrent, 
  PollardRhoBrent/ext, QuadraticSieve, SmallFactors, 
  SmallFactors/ext, SmallFactors/lib, anm1, easy, easyfunc, 
  from_signature, ifact0th, ifact1st, ifact235, lenstra, mixed, 
  morrbril, mpqs, mpqsmixed, noext, noext/ifact0th, 
  noext/ifact1st, pollard, pollp1, power, pp100000, squfof, 
  wheelfact]

Now we can test them all. A preliminary version of the following test indicated to me that numbers 14, 15, 16, 24, 25, and 29 were difficult to test: They either produced unrecoverable errors on my test case, or they went on way too long. So, here's the test code:

n:= nextprime(10^24)*nextprime(10^25):
for k,M in Meths do
    if k in {14,15,16,24,25,29} then next fi;
    try
        print(k, M);
        print(CodeTools:-Usage(ifactor(n, M)))
    catch:
    end try
od:

From doing this, I couldn't find any evidence of multiple-core usage. That doesn't mean that it's not there, only that I didn't see it. Anyway, the above should give you a framework for further exploration.

A brief tutorial on "display-form" inertness in Maple

Maple offers several forms of inertness, i.e., ways to control when, if ever, expressions or particular subexpressions evaluate. I'm not aware of names for these different forms of inertness. So, for now at least, I'm going to call the topic of this thread and this article "display-form" inertness because it's primarily (although not exclusively) used to manipuate the prettyprinted displayed form of expressions while keeping them within manifestly mathematical data structures that--unlike for example strings--are only slight variations of their noninert (or active) counterparts. This form of inertness is mostly controlled by using the command value, by prefixing symbols and operators with %, and by the package InertForm. The package Typesetting and procedures whose names begin `print/ are also involved, although this aspect is usually transparent to the end user and often also transparent to the programmer who's one step removed from the end user.

The option inert=true vs. inert=false in package InertForm: The inert=true and inert=false options to some commands in InertForm do not change whether the expression is inert; they only change whether the inert operators are displayed as inert (i.e., in gray) rather than as normal (i.e., in blue).^[*1] And it's only the inert operators, not the entire expression, whose color changes.

So now that I've pointed out the subtlety of this difference, do you notice it? It's clearly visible in the above Answers by @acer and @Kitonum .

If neither of these options is used, that's equivalent to inert=true.

Making symbols inert with %: Any Maple symbol (e.g., command name (if it's a procedure), function name, variable, named constant (e.g., Pi, infinity)) can be made inert by prefixing an %. Although these inert forms usually "do nothing" (i.e., they are truly "inert" as in the common usage of that word), that is not actually a requirement. Executable code in procedures can be supplied for them, and those procedures can return an inert form as an unevaluated function. One possible use for that is to do syntax checking of arguments without doing the computation. An example of this is shown in the Reply below.

Special-character quoting (e.g., `%Pi`) is usually not needed (in 1D input (I don't know about 2D Input)) when prefixing %, but it is allowed. Multiple levels of inertness can be applied by prefixing multiple %. In that case, the quoting is required.

Inert arithmetic infix operators: The five basic^[*2] infix arithmetic operators as well as noncommutative multiplication are available in inert forms for infix use: %^, %., %*, %/, %+, %-. These can be used infix in 1D input. I don't know about 2D Input, but your experience mentioned above suggests that they can't be. As Kitonum showed, they can always be used in prefix (aka functional) form, in which case they're just symbols, and the previous paragraph applies. You cannot create multiple levels of inertness for the infix forms.

You cannot use quotes with any infix operators (regardless of whether they're inert) if they're being used in infix form. The quotes ` ` (sometimes called "back quotes") turn anything into a symbol. Infix operators cannot be symbols, although almost all of them have an equivalent symbol form that allows them to be used in functional form (just like Kitonum used `%*` above).

Like any "multiple-punctuation-mark" infix operator in almost any computer language^[*3], you cannot separate the characters with any white space. And of course ab or `ab` (those are identical) and `a b` are distinct symbols.

Other historically inert commands: The commands Eval, Diff, Int, Intat, Limit, Sum, Product, Normal, and Hypergeom (and perhaps I've missed some) are inert forms of the corresponding uncapitalized commands. These are much older than is using % for inertness; however, their form of inertness is essentially the same as with %.

There are many other cases of two Maple commands whose names only differ by capitalization and the capitalized form is inert; yet the lowercase form is not an active form for it. In these cases, the capitalized form is usually a data container for input in some specialized domain of computation such as mod, modp1, and evala. See their help pages for examples of these inert functions.

Removing inertness with value: The value command applied to an expression removes one level of inertness from every inert symbol (including the historical ones from the previous paragraph) and inert infix operator in the expression. If special techniques need to be used to remove the inertness, value has a facility for supplying procedures^[*4] for that. Likewise, any otherwise unused symbol, regardless of any special characters, can be made into an inert form of any other symbol (i.e., there need not be any connection between the spellings of the inert and active forms) by teaching^[*4] value about it. The previous sentence applies to symbols only, not to infix operators. It is not necessary that any noninert form produced by value actually does anything, i.e., it could also be "inert" for all practical purposes.

The InertForm package: In summary, there are numeous ways to use inertness without using the InertForm package. However, one way (among several) that the package is quite useful is for working with expressions that are so volatile that they must be entered as strings before being Parsed by the package into mathematical expressions. Examples of such volatile expressions are those that would otherwise be changed by Maple's numerous forms of automatic simplification.

The command Typeset produces displayable expressions whose detailed internal structure can be examined with lprint. This is useful for learning how the Typesetting package is used (essentially as a mark-up language).

`print/...` procedures: A procedure named `print/...` where ... is a function name (usually an essentially inert function, although its inertness need not be formally imposed by % or by any other means) is one of the easier ways to create customized prettyprinted displays for inert functions. For example, the inert function can be broken down into simpler inert functions whose prettyprinted display is already handled by the system. An example of this is in the Reply below. 

These procedures could manipulate the output of InertForm:-Typeset. An example of this---performing a small "surgery" to correct the blue exclamation points---is also in the Reply below.

Relevant help pages: The most-important help page for more details on what I've presented here is ?value. Also relevant are ?operator,precedence, ?print, and of course ?InertForm (in particular, ?InertForm,Display and ?InertForm,Typeset). Some more details about using "punctuation mark" operators as function symbols can be found at ?use.

The commands FromInert and ToInert are not related to the "display-form" inertness being discussed in this thread and this article. They are oriented towards debugging, automatic code generation, and what I call "code surgery".

I've used several technical terms in this article without providing definitions. These appear in upright bold if they're available verbatim as Maple commands, keywords, etc., or in italics if they're being used in a technical sense (such as specific to formal computer science) rather than as common English. Feel free to ask me for specific definitions of any of these.
 

[*1] These styles, such as the colors gray and blue, could be changed if you really wanted. I've never seen it done, and I'd recommend against it unless you were a master craftsman at typography and document layout, lest you significantly reduce the readability of the output, especially its readability to experienced Maple users who are expecting certain styles.

An example in the Reply below corrects a bug that causes the color of the exclamation points to be default blue rather than inert gray in inert factorials.

[*2] Personally, I'd also consider mod to be a "basic" infix arithmetic operator, but it doesn't have an infix inert form. Just like most of the infix "punctuation-mark" operators, its inert infix effect can be achieved by quoting, i.e., by using `%mod` as a function symbol.

[*3] Fortran, at least in some of its older forms, is unlike most languages in that it has some exceptions to these white space restrictions.

[*4] Specifically, this is done by defining a procedure named `value/...` where the ... is the inert function symbol. The details are given at ?value. An example is shown in the Reply immediately below.

There's an option to dsolve whose sole purpose is to get around this specific error. It's called continuation. See if you can find the help for it. It's probably on the help page ?dsolve,bvp_numeric,advanced (that's off the top of my head, so, not sure about that page name). If you haven't figured it out in an hour or so, I'll post an Answer. 

Also, I've answered this exact question on MaplePrimes many times. If you're lucky, you may be able to find such an Answer. Unfortunately, I can't, because the Search Tool here is horrible. It would likely take several days to hand scan through all my Answers looking for the ones related to numeric solution of boundary-layer BVPs.

Good Question; Vote Up.

There are a few snags in your formulation. These are not quite mathematical issues but rather issues of the level of symbolic abstraction, symbolic simplication, and when (if ever) variables will be given numeric values. My solution below assumes that n be specified up front (as you have done), and i be specified before the merger of the piecewises. The result for each i is a single piecewise function that is branched wrt x/h not x. So, it doesn't quite make sense that phi[i] be specified as a function of x, because this process requires x and h to be names while phi is operating. Note that the result can be expressed as a function of x/h, which for a piecewise makes it much simpler than a function of x and h considered independently. The result of phi[i](x) can be applied to numeric values of x and h. I chose to do it this way so that the final result would be most amenable to further symbolic processing, such as integration. If instead the goal was a crude numeric application such as plotting, these subtler details wouldn't matter much. Let me know if you'd like any of these details rearranged.

Download PiecewiseMerge.mw

I can't post the worksheet inline, but you can download it. Here is the code portion:

restart
:
S:= x-> piecewise(
     x <= -2, 0,
     x <= -1, (2+x)^3/4,
     x <= 0,  (2+x)^3/4 - (1+x)^3,
     x <= 1,  (2-x)^3/4 - (1-x)^3,
     x <= 2,  (2-x)^3/4,
     0
):
n:= 9:
_phi:= proc(x::name)
local 
    i, xi, h, R,
    T:= table([
        0=    S(xi) - 4*S(xi+1),
        1=    S(xi-1) - S(xi+1),
        _n=   S(xi-_n) - S(xi-(_n+2)),
        _n+1= S(xi-(_n+1)) - 4*S(xi-(_n+2))
    ])
;
    if not procname::indexed or nops(procname) <> 1 then
        error "phi must be called with a single index"
    fi;
    i:= op(procname);
    if not i::integer[0.._n+1] then
        error 
            "phi's index must be an integer from 0 to %1 but"
            " received %2", _n+1, i
    fi;
    R:= `if`(assigned(T[i]), T[i], S(xi-i));    
    unapply(subs(xi= x/h, convert(R, piecewise, xi)), [x,h])
end proc
:
phi:= subs(_n= n, eval(_phi))
:
#Examples:
for k in [0,1,iquo(n,2),n,n+1] do print(i=k, phi[k](x)(x,h)) od:

​​​​​

Maple 2021 has introduced the reduce option to seq. Replacing add with seq[reduce= `+`] is very efficient, and it avoids the kernel crash caused by add. In the code below, I have made your B a procedure, like Tom did. But using the reduce option means we don't need a list or to use elementwise operation (B~). This is thus more efficient (especially memorywise). (And add would be efficient also, but there's a bug at the moment.)

I've simplified and generalized your expression B such that A, B, and any higher order such expression can be handled by the new procedure B below. As you can see, the time to evaluate B over all tuples and add them is barely measurable.

restart:
Iterator:-CartesianProduct([1,1]): #Force compilation

P:= 2..4: #suffixes of m and c variables
N:= [$0..4]: #evaluation values of m variables

B:= subs(
    {_C= [$P], _V= [c||P]}, 
    proc(M)
    local r:= add(M*~_C), s:= 1+r, t:= s-add(M);
        r!*c0^t*mul(_V^~M)/t!/mul(M!~)/c1^s
    end proc
):
CodeTools:-Usage(
    seq[reduce= `+`](B(v), v= Iterator:-CartesianProduct(N$(rhs(P)-1)))
):
memory used=0.69MiB, alloc change=0 bytes, 
cpu time=0ns, real time=3.00ms, gc time=0ns

Tables are one of the fundamental data structures built directly into the Maple kernel. They are extremely efficient and extremely flexible. They are quite suitable as the foundation for building numerous other data structures. Indeed, they could easily be used to implement MultiSets (and perhaps they are so used).

Is this suitable?:

CodeGeneration:-Matlab(
    proc(a::Vector, N::posint)::hfloat; 
    local i::posint; 
        Sum(a[i], i= 1..N) 
    end proc,
    resultname= "sum"
);
function sumreturn = sum(a, N)
  r = 0;
  for i = 1:N
    r = r + a(i);
  end
  sumreturn = r;

All the ::... declarations could be omitted. An arrow procedure could be used:

CodeGeneration:-Matlab(
    (a, N)-> Sum(a[i], i= 1..N), 
    resultname= "sum"
);
Warning, procedure/module options ignored
function sumreturn = sum(a, N)
  r = 0;
  for i = 1:N
    r = r + a(i);
  end
  sumreturn = r;

The warning can be ignored in this case.

What you're trying to do is trivial with Iterator:-CartesianProduct:

P:= 4: M:= [m||(2..P)]: N:= 4:
add([
    for c in Iterator:-CartesianProduct([$0..N]$(P-1)) do 
        eval(B, M=~[seq](c)) 
    od
]);

sum2N:= proc(N::nonnegint) 
local summation:= 0, i; 
    for i from 0 to N do summation:= summation+(N+i)^2 end do
end proc;

The technique that I'm about to provide will work for any numeric dsolve solution, regardless of the number of equations, the order of the derivatives, or whether it's IVP, BVP, or DAE. First, the ODEs need a slight rearrangement. In your case, that's

deq1:= diff(y(t), t, t)= -y(t)*abs(y(t));

All that I've done is solve the ODE for its highest-order derivative. Once that is done, any derivative of y of any integer order j (including j=0 or j=1), can be expressed as 

eval[recurse](diff(y(t), [t$j]), {deq1})

and these expressions will be recognized by odeplot. The only modification to the code that is needed is to solve the ODEs for their highest-order derivatives. So, for example,  you can make an array plot of y and its first 4 derivatives by

deq1:= diff(y(t), t, t)= -y(t)*abs(y(t)); 
ic1:= y(0) = 1.0, D(y)(0) = 0.; 
dsol1:= dsolve({deq1, ic1}, numeric, range= 0 .. 10); 

plots:-display(
    Array(
        0..4,
        j-> plots:-odeplot(
            dsol1,
            [t, eval[recurse](diff(y(t), [t$j]), {deq1})],
            labels= [t, ``],
            title= typeset(diff(y(t), [t$j]))
        )
    )
);

In your first piecewise command, you've misspelled eps as esp. In your second piecewise, you've done that again, and you've misspelled beta as betta.

In  your last command, t_atan2 means nothing to Maple, unless you've defined it elsewhere. You should replace it with arctan. Maple's arctan can process 1 or 2 arguments.

Using the command subs as you've shown doesn't "set" anything; it can only do so if its final argument is some expression that contains beta.

I recommend against using and (or other logic operators) in the conditions of piecewise. The processing of the conditions proceeds left to right, thus each subsequent condition automatically includes the negation of the disjunction of all previous conditions. Making use of this fact, you almost never need to use and, etc. In other words, 
piecewise(cond1, expr1, cond2, expr2, cond3, expr3)
is logically equivalent to
piecewise(
    cond1,  expr1,
    not cond1 and cond2,  expr2,
    not (cond1 or cond2) and cond3,  expr3,
    not (cond1 or cond2 or cond3),  0
);

So, I recommend this:
W:= beta-> piecewise(
        beta < Pi/2-theta__M, 
            0,
        beta <= Pi/2+theta__M, 
            (1+epsilon*cos(4*theta__M))*cos(4*beta)/cos(theta__M),
        beta <= Pi-theta__M, 
            1+epsilon*cos(4*beta)
);
W(arctan(diff(pho(x,y),y), diff(pho(x,y),x)));
 

What you're calling a "comma-separated vector" Maple simply calls a list. Your x and f are lists. There is no need for you to use Matrix or Vector for what you show above; it can all be done with lists:

x:= [x__1, x__2, x__3]:
f:= x__1: #no need for f to be a list
Gradf:= diff~(f, x);
Hessf:= map2(diff~, Gradf, x); 

What you call Gradf2 is usually called the Hessian of f.

If you still want to use matrices or vectors, a matrix or vector v can be converted to a list by [seq](v).

Don't use evalf. Its sole purpose is to make decimal approximations. Here are a few ways to do what you want:

1. exp(2*Pi*I*i/3) $ i= 0..2
2. seq(exp(2*Pi*I*i/3), i= 0..2)
3. exp~(2*Pi*I/3*~[$0..2])

As we have discussed here several times in the past, you are not using mod correctly. Use irem instead of mod.

You have found a bug in Maple, and I'll give you a workaround in a moment. But your code itself also has 1 serious bug and 1 serious shortcoming, and it amazes me that these didn't confound the issue of the true bug that you found.

Bug: Change
f:= x-> piecewise(0 <= t, 1, t < 0, t - 1)
to
f:= t-> piecewise(0 <= t, 1, t < 0, t - 1)

Shortcoming: Change
F:= x-> int(f(t), t= 1..x)
to
F:= unapply(int(f(t), t= 1..x), [x])

Workaround: Change the unapply to
F:= unapply(int(f(t), t= 1..x, method= FTOC), [x])

So, why use unapply? Because F contains a symbolic operation, the integration, that can be performed independently (as well as once-and-for-all) of F's parameter x.