@emendes 

Since the vast majority of the cputime used for this computation by all the methods shown so far is for the garbage collection, and given that Maple's garbage collection is already parallelized, it might be worthwhile to increase that parallelization. Please redo with SubsetPairsSparse the timing of one of the large computations that you've already shown with this one change: Before starting the timing, give this one command:

kernelopts(gcmaxthreads= kernelopts(numcpus));

@vv Maple has a massive and very efficient infrastructure for doing set problems such as this one. To write this efficiently in a compiled language such as C would require far more effort than "reinventing the wheel". Note that for all the large computations posted here, the majority of the time is in garbage collection, which is automatically parallelized in Maple. 

@Preben Alsholm 

Vote up.

I'd like to point out to the OP that the mathematical and computational essence of what you just said is the same as what I said. Of course, you've also given a detailed example. 

@Robert Matlock 

I apologize that the rhetoric of my comment was flippant; however, I stand by its mathematical and computational content 100%.

I have read and fully understood what you wrote as well as all previous MaplePrimes posts on this subject, and I've used those techniques many times, both for IVPs and BVPs. My comments here apply only to the case of evaluating the highest-order derivatives that actually appear in the system, not to other derivatives or integrals. There does not exist a more-accurate way (i.e., a way with more error control) to numerically evaluate those derivatives than direct numeric evaluation of the necessarily-already-present symbolic expressions that specify those derivatives. The error control built in to the numerical arithmetic evaluators (evalf, etc.) is already necessarily at a higher precision than that used by dsolve.

It was the truth (and, to me at least, the obviousness) of the "necessarily" in that last sentence (as well as a toothache and lack of sleep) that led to the flippancy of my previous Reply.

@Robert Matlock Suppose that dsolve did provide the highest-order derivative in its output. By what magic process would it compute that other than the way that we're describing to you? 

@emendes 

Please clarify the problem's size for the large example that you just ran. You referred to the "length" of the two input lists (which I'll call A and B throughout this thread). If this refers to Maple's length command, then that information is not very useful. Knowing the problem's "size" (via multiple metrics) is essential for guessing which algorithm will be most efficient.  So please execute the following, which should run quickly even for very large input.

ProbSize:= (A::list(set), B::list(set))->
    [nops~, (add@nops~)~, (nops@`union`@op)~]([A,B])
:
#For example, using your initially presented problem:
ProbSize(ans6, ans7); 
             [[20, 30], [120, 210], [16, 14]]

And if you still have access to the output of your test run (which I'll call L), please do the following:

S:= (add@nops~)(L); 
            S:= 47

The number S provides a theoretical lower bound on the time complexity of the computation.

Naive algorithms for this problem (such as yours, vv's, Tom's first, and mine above) have a runtime proportional to nops(A)*nops(B). (The storage required by vv's (but not the others) is also proportional to this, unfortunately.) For the example that you initially presented, this product is 600; for the example you just tested, it's about 6 billion (yikes!) if by "length" you meant nops.

To not be naive, let's consider the sparsity of the problem. The numbers that I just asked you to compute are attempts to measure the sparsity. Comparing 47 to 600, the example problem is quite sparse. Tom's second algorithm attempts to exploit the sparsity (whether he knew it or not). But his Maple implementation of the algorithm has a serious flaw that makes that procedure nearly worthless for large input. The flaw is only in the implementation, not the algorithm itself. That flaw is the loop 

for i in L1[j] do T[i]:=`union`(T[i],{j}) od

As is well known and often discussed here, it's extremely inefficient to construct a long list, sequence, or set by iteratively appending, the time to do so being proportional to the ​​​​​​square of the number of iterations.

The problem in Tom's code can be fixed, and here is my version that does that:

SubsetPairsSparse:= proc(A::list(set), B::list(set))
local 
    a, b, x, j, 
    T:= table([for x in `union`(A[],B[]) do x= table() od])
;
    for j,b in B do for x in b do T[x][j]:= () od od;
    T:= {indices}~(T, 'nolist');
    [for a in A do `intersect`(seq(T[x], x= a)) od]
end proc
:

Notice that neither this nor Tom's version of it do any looping process (for, seq, map, ~) whose number of steps is nops(A)*nops(B). Please try running that on a largish example, compare with my SubsetPairs1a, and report back. If the example that you use is not the one mentioned in your most-recent Reply, then please run my problem-sizie metrics on it also. If that goes well, we'll work on a threaded version.

Regarding the bipartite graph: Sorry, I forgot to relabel the vertices, which is needed for it to be recognized as bipartite. I've corrected the code above, and the entries from A will appear as negative numbers on their vertex labels. I also added to the DrawGraph command some basic options to improve readibility. Much more could be done along these lines.

Given an instantiation to numeric values as shown above, you may want to automatically verify whether those values are all real. That can be done like this (starting where the Answer above left off):

Inst:= eval(BoundV, FreeV) union FreeV;
type(Inst, set(name=realcons));
                              true

Can you provide a web page or paper that gives details of the algorithm? The information that I can find so far is too general to be of much use.

Just to clarify and to make sure: You're interested in an algorithm to factor multivariate polynomials with only rational coefficients, right?

Here is a complete Maple implementation covering all documented calling sequences of the Matlab function, and allowing Y to have any number (including 0) of dimensions, singleton dimensions, and starting indices other than 1. Of course, the vast majority of the code below is to handle all the weird cases; the mathematical computation of the discrete trapezoid rule is just a dot product.

trapz:= proc(x::rtable, y::{rtable, nonnegint}, dim::nonnegint, $)
description "Maple implementation of Matlab's trapz";
option `Author: Carl Love <carl.j.love@gmail.com> 2020-Dec-9`;
uses AT= ArrayTools;
local X, Y, Dim, nd, s, dims;
    #Deal with the 1-, 2-, or 3-argument calling sequences, and decide between
    #the 2 possible cases for the 2-argument sequence.
    if nargs=1 then Y:= x
    elif nargs=2 then 
        if y::rtable then (X,Y):= (x,y) else (Y,Dim):= (x,y) fi
    else (X,Y,Dim):= (x,y,dim)
    fi;

    #Get Y's dimensions:
    if (nd:= nops((dims:= [rtable_dims](Y))))=0 then return fi;
    local sz:= (rhs-lhs+1)~(dims);
    local LB:= (1-lhs)~(dims);

    #Find the key dimension of Y if it wasn't specified in the
    #arguments:
    if Dim::name then 
        if Y::Matrix then (Dim,s):= (2, sz[2])
        #Find 1st non-singleton dimension:
        else for Dim,s in sz while s=1 do od
        fi
    else
        try s:= sz[Dim]
        catch: error "Y has no non-singleton dimension"
        end try
    fi;

    #If X was passed in, check its size, and adjust its shape if
    #necessary; if not, set it to be unit spacing:
    if X::name then X:= `<|>`(1 $ (s-1))
    elif rtable_num_dims(X)<>1 or numelems(X)<>s then
        error "expecting X to be 1-D with %1 elements", s
    else 
        X:= AT:-Alias(X, [s], ':-row');
        X:= X[2..] - X[..-2]
    fi;

    #Prepare Y to be handled essentially as an Array of column vectors
    #(without actually breaking it up):
    local Yd:= AT:-Alias(Y/~2, sz, if nd=1 then ':-column' else fi);
    Yd:= 
        Yd[(()..()) $ (Dim-1), 2.., (()..()) $ (nd - Dim)] +~
        Yd[(()..()) $ (Dim-1), ..-2, (()..()) $ (nd - Dim)];
    s--;
    LB:= LB[[1..Dim-1, Dim+1..-1]];

    #final computation:
    if nd=1 then
        X.Yd
    else
        rtable(
            dims[[1..Dim-1, Dim+1..-1]][],
            (()-> 
                local J:= [args] +~ LB;
                X . AT:-Alias(
                    Yd[J[..Dim-1][], .., J[Dim..][]], [s], ':-column'
                )
            ),
            ':-subtype'= 
                if Y::Matrix then Vector[['-column', ':-row'][Dim]] 
                else Array
                fi
        )
    fi
end proc
:

@acer Okay, I understand your point. I was thrown off by the OP's use of the word precedence. Now that I reread the Question in light of your Reply, I see that the issue has nothing to do with precedence, and that the OP had simply misused that word (albeit in good faith).

@wswain You wrote:

• The original array elements were of a text(string) type, ...

Yes, that's correct.

• ... but you use subsindets to "transform" the list elements to string. 

No, that's incorrect. The second argument to subsindets, string in this case, tells the command to look for objects of type string in its first argument, the list. Hence, the list must already contain strings for this to work.

• Can you explain the "transformer" and "rest" arguments ...

The "rest" argument is not used in this example. If it had been used, it would be a constant (often an option) that would be the fourth argument of subsindets and passed as the second argument to the transformer.

• ... how you know a valid transformer operator (in this case 'string') ,,,

You mean operand, not "operator". The subsindets command only passes to the transformer objects of the correct type, string in this case, specified in its second argument.

• ...  vs it just acting on the elements?

Indeed, subsindets is overkill in this flat, linear case. It could be replaced (in Maple 2018 or later I believe) by

map(index, [seq](Array3[1, 2..]), [5..])

•  I presume a different substring operation could have been used to remove chars 1-4 ?

No, u-> u[5..] is precisely the substring operation that does remove characters 1-4.

The Answers by Tom and acer show that some relatively minor bug related to your specific Question was fixed by Maple 2016. However, the more-important bug, the operator precedence itself, was not fixed until Maple 2020 (I think at my request). Here's proof:

kernelopts(version);
  Maple 2019.2, X86 64 WINDOWS, Nov 26 2019, Build ID 1435526
a %* b %+ c %* d;
                     %*(%+(%*(a, b), c), d)

kernelopts(version);
  Maple 2020.1, X86 64 WINDOWS, Jul 30 2020, Build ID 1482634
a %* b %+ c %* d;
                     %+(%*(a, b), %*(c, d))

@Madhukesh J K At the moment, I think that you'll learn more by trying to follow written directions rather than an explicit code example. If you try this, and it doesn't work, please post an executed worksheet showing your attempt, and I'll give you a code example.

First, do you know what a coefficient is? And do you know what a nonlinear term in a differential equation is? If you're not sure about either of these, please ask. Otherwise, here are the explicit and simple steps:

1. Make C a coefficient of every nonlinear term in all 4 ODEs. In other words, multiply each such term by C.
2. Add option continuation= C to the dsolve command. 

That's all there is to it.

@Wilks You have taken a very easy problem and made it enormously complex. Like I said from the start, there's no need to use solve. Like I said after you provided details, there's no need to use int or intat. All that can be replaced by dsolve.

You need to put restart at the beginning of your worksheet.

The entire "chapter" of your worksheet entitled INTEGRACIÓ EQ. DE L'ELÀSTICA (LLEI DE DESPLAÇMENTS) can be replaced by these few lines of code, which computes all 8 Ys at once:

for Y in {Y||(1..8)} do
    unassign(Y); #just in case
    assign(
        Y= unapply(
            (rhs@dsolve)
                ({diff(Y(x), x$2) = cat(Y,`"`)(x), Y(0)=0, D(Y)(0)=0}),
            x
        )
    )
od:

That's if you want the Ys without the symbolic constants of integration. If you do want the constants, it can still be done with dsolve; let me know. It seems like you want the same Y's to have the constants in one part of the worksheet but not the other. Just use different names to distinguish them!

@David Sycamore Given any list L of [p,q] pairs (such as is generated by Kitonum's or my procedure), the sequence of just p can be extracted as Pseq:= op~(1, L)[]. The efficiency that can be gained by not storing the q in the first place is minimal.