This is the appropriate place to discuss MaplePrimes problems. That's why there's a "MaplePrimes" entry in the Question-header "Product" check boxes. 

Additional details: I have two computers---both Windows 10, but different versions---and I have this problem on both. For both, the browser I use for posting is Google Chrome.

@nm If your guess is correct, then I've been extraordinarily unlucky. Like I said, I haven't been able to display any worksheet at all for about a year although I've tried many times.

@Ronan The approach that you show is a last resort for infinite sums for which evalf won't work. Vv's Answer shows that that isn't the case here.

Vote up. But due to the presence of LinearAlgebra in your title, I'd like to clarify something that may not be obvious to all readers: Your examples display the same (whether via Latex or GUI) regardless of whether they are created with LinearAlgebra. Indeed, many users have the false impression that LinearAlgebra is required to use Vector or Matrix. No, those are top-level commands.

I think that one of the examples that you present should be condidered a bug, just for the sake of consistency.

I dispute your contention that the n returned by boo is local to boo because of this:

restart;
foo:=module()
   export n::integer:= 2;
   export boo:=proc()::integer;
     return n;
   end proc;       
end module
:
foo:-boo();
                               2

Also consider this:

restart:
foo:= module() export n; end module:
type(foo:-n, `local`);                
                              true

So, it is the presence of the ::integer that makes the difference.

@tomleslie The last algorithm that you show impressively exploits the problem's sparsity. There was an efficiency problem with your Maple implementation, for which I've given a correction in my Answer. Nonetheless, I wanted to directly congratulate you here for coming up with a great algorithm, which is the most difficult creative aspect.

@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
: