Here's my procedure for the lexicographic product of graphs. I made it work for any number of factors.

LexProd:= proc(G::seq(Graph), $)
uses GT= GraphTheory;
local
    LexProd2:= proc(G::Graph, H::Graph)
    local 
        (Vg,Vh):= op(GT:-Vertices~([G,H])),
        (ng,nh):= op(`..`~(1, nops~([Vg,Vh]))),
        (Ng,Nh):= op(op~(4, [G,H])),
        i, j, J,
        P:= [seq](seq([i,j], j= nh), i= ng),
        k:= 0, K:= table((p-> op(p)= ++k)~(P))
    ;
        GT:-Graph(
            (curry(nprintf, "%a:%a")@((i,j)-> (Vg[i],Vh[j]))@op)~(P),                
            (
                ((i,j)-> {
                    seq(seq(K[k,J], k= Ng[i]), J= nh), 
                    seq(K[i,k], k= Nh[j])
                })
                @op
            )~(P)            
        )
    end proc
;
    if nargs=0 then error "at least 1 graph needed"
    elif nargs=1 then G[1] 
    else foldl(LexProd2, args) 
    fi
end proc
: 

If you cube both sides of the equation to remove fractional exponents, then you'll get 3 solutions, 2 of which are trivial, and the 3rd is much simpler than your solution above. Using odetest~ quickly returns 0 for all 3.

restart;
ode:= diff(y(x),x)^3 = x^3*(y(x)^2-1)^2;
sol:= [dsolve](ode);
odetest~(sol);

If you just want to extract column j as a vector from matrix M, you could use M[.., j]. And unlike the Column command, this simple indexing also works on the left side of the assignment operator.

You are correct that the underlying issue is closely related to MutableSets being object modules. The "recommended way" would be for Maplesoft to include a ModuleType procedure in the object definition of MutableSet. It's a shame that they haven't, and it seems ridiculous to not do so for any container object (such as MutableSet).

You can use this as a workaround:

M:= MutableSet(a, b, c);
type(M, And(MutableSet, set(symbol) &under (convert, set)));

In a procedure header, that could be

proc(M::And(MutableSet, set(symbol) &under (convert, set)))

In either of the above, replace symbol with your desired type for the elements. 
 

Your method of writing to a file which you then reread to get the result of a simple computation that's been displayed as userinfo is extremely convoluted, slow, and low precision (4 digits). Instead, just do that simple computation, which in this case can be done by the command LinearAlgebra:-NullSpace. So, you have a 3-column matrix named blu. You consider each row of blu to be the normal vector of a plane. You want a basis for the vectors in that plane. If you want to do this for the i^th row of blu, then the set of 2 basis vectors is simply

LinearAlgebra:-NullSpace(blu[i]);

This can be done for the entire matrix blu with a single command. There are numerous ways to do that depending on the exact output format that you want. For example, the following returns two 3x200 matrices B1 and B2 such that the i^th column of B1 is the first basis vector corresponding to blu[i], and that of B2 is the second basis vector:

(B1,B2):= map(`<|>`@op@op~, [1,2], LinearAlgebra:-NullSpace~(convert(blu,listlist)))[];

That command uses 0.12 seconds on my computer to do all 200 rows.

There are a few reasons why it doesn't work, and several ways to fix it. I'm guessing that you intend for x1, x2, x3 to be the roots of eq? But you never substituted them in for x in eq. I think this is the simplest fix that does what I think you really want:

restart:
eq := x-> a*x^3 + b*x^2 + c*x + d;
 s := x1 + x2 + x3 = -b/a;
 sp := x1*x2 + x1*x3 + x2*x3 = c/a;
 p := x1*x2*x3 = -d/a;
sol:= ()-> solve({eq~([x1,x2,x3])[], p, s, sp}, {x1, x2, x3},explicit);
sol(); #no solution, yet                       
 a := 1;
b := -5;
c := 6;
d := -1;
 sol(); #very long solution quickly returned   

It seems to be a bug such that the columns of the selector matrix A are used in reverse order. You can use this instead:

eval~(`if`~(A, B, 'NULL'));

If I ignore the typeset formulas in your Question---which necessarily have errors as alreadyb pointed out---and I focus instead on the text of your Question, then I think that the command BesselJZeros is the major part of the solution. See help page ?BesselZeros.

The number of terms isn't a mathematically invariant property of expressions. In other words, it's not a property that's necessarily preserved by transformations under which expressions remain mathematically equal. Rather, it's a structural property. So, expand shouldn't be used. This is sufficient:

nTerms:= (e::algebraic)-> `if`(e::`+`, nops(e), 1)

Here's a procedure that returns the values of the four formulas in your Question. This code likely requires 1D Input. I can modify it if need be, but I strongly recommend learning 1D Input.

ANOVA_3D:= (X::And(rtable, 3 &under rtable_num_dims))->
local
    N:= numelems(X), 
    Add:= ArrayTools:-AddAlongDimension,
    Xij:= Add(X, 3),
    SS:= A-> add(A^~2)*numelems(A)/N,
    SSi, SSj, all    
;
    Record(
        "SS__P"=  (SSi:= SS(Add(Xij, 2))) - (all:= SS(<add(Xij)>)),
        "SS__A"=  (SSj:= SS(Add(Xij, 1))) - all,
        "TSS"=    SS(X) - all,
        "SS__AP"= SS(Xij) - SSi - SSj + all
    )
:
#Example usage:
(n, k, r):= (90, 100, 110): #There'll be nearly a million entries.
#essentially the same as Array, but allows the random initializer:
X:= rtable((1..n, 1..k, 1..r), frandom(0..2, 1), datatype= hfloat):
R:= ANOVA_3D(X);
                                                               
R := Record(SS__P = 35.6539282805752, SS__A = 35.7316001050640, 
  TSS = 330279.334091332, SS__AP = 2965.63736086327)

#Access values like this:
R:-SS__P;
                        35.6539282805752

Let me know when you want to access values from the F-distribution.

You'll need 3 initial conditions. I used 3 random values in interval 0..1. This system of ODEs is known to be particularly sensitive to changes in initial conditions; the term "butterfly effect" was essentially invented because of this system.

For the plot, I used thickness= 0.1. A very fine line is ideal for this. But you'll need Maple 2022 for that. If you have earlier Maple, change 0.1 to 0. Maple has numerous schemes to color the plot.

restart:
Digits:= 15:
V:= [x,y,z]:
eqs:= diff~(V(t), t)=~ [35*(y-x), -x*z-7*x+28*y, x*y-3*z](t);
sol:= dsolve({eqs[], (V(0)=~ V||~0)[]}, numeric, maxfun= -1, parameters= V||~0);
sol(parameters= ['rand(0.0..1.)()' $ 3]);
plots:-odeplot(sol, V(t), t= 0..50, numpoints= 5000, thickness= 0.1);

To get the first entry from each pair (where "first" is with respect to the order displayed, which is not necessarily the order that they were entered) do

op~(1, ({entries}@op~@ListTools:-Classify)(lhs, AA));

This should work in Maple 12:

c:= plots:-spacecurve([cos(t), sin(t), t], t= -Pi..Pi);
map2(op, 1, indets(c, specfunc(anything, CURVES)))[1];
 

Like this:

restart:
oldGAMMA:= eval(GAMMA):
unprotect(GAMMA):
GAMMA:= overload([
    proc(x::And(integer, Not(0)), y::And(algebraic, Not(complexcons)), $)
    option overload;
        'procname'(x,y)
    end proc,
    oldGAMMA
]):
protect(GAMMA, oldGAMMA)
:
#Test:
GAMMA(-1,x), oldGAMMA(-1,x), GAMMA(3);

Like this (I've divided all the quantities by 1000 to reduce visual clutter):
 

Download LPTransNetw.mw