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

@Earl (Perhaps that title should be "A get-Maple-to-do-more-algebra approach.") I should've shown that simplification of W. I think that that would've avoided that version difference of implicitplot. And it'll make your attempt to solve W for y work better. Indeed, by solving for x and y separately and plotting both, you can get a perfect plot (shown below).

We simplify W by using the fact that it's a rational function of x and y implicitly equated to 0. Thus,

1. Its denominator doesn't matter and can be discarded. This is the purpose of numer below.
2. The "content" (common factors of all numerator terms) can be discarded. This is the purpose of op&<2 in 
   op~&<1 @ op&<2 @ evala@AFactors below.  (My &< is defined below.) That @-expression is equivalent to
   e-> op~(1, evala(AFactors(e))[2]).
3. Exponents on factors of the numerator can be discarded. This is the purpose of 
   op~&<1. 

restart
:
#abbreviations for things that I commonly use in @-expressions:
`&<`:= curry: `&>`:= rcurry:   
:
WattEq:= r^2 - (b^2 - (a*sin(theta) + sqrt(c^2 - a^2*cos(theta)^2))^2);


#Don't use decimal values for the parameters.
params:= (`=`~ &< [a,b,c])~([[31,11,30]/10, [1,sqrt(2),1], [21,22,6]/10]);


rng:= -b..b:
V:= [x,y]:
abc:= e-> eval(e, params[Params]):
 
Params:= 1: #which parameter values to use

W:= (op~&<1 @ op&<2 @ evala@AFactors @ numer @ evala@Norm @ eval[recurse])(
    abc(WattEq),
    #polar to cartesian: 
    [r= sqrt(x^2+y^2), (cos,sin)(theta)=~ V[]/~r]
):
if nops(W)<>1 then WARNING("multiple factors") fi:
W:= mul(W);


plots:-implicitplot(abc([W, V[]=~ rng])[], scaling= constrained);

For the purely functional/algebraic plot, we continue with the same W, and solve it for x and y separately. This'll lead to much duplication of the plotted points, but that doesn't matter; the gaps left by each will be filled by the other. The solutions for x are plotted in parametric form (e.g., plot([f(y), y, y= a..b])) so that the axes are switched. The solutions are substituted into a procedure template that does hardware-float evaluation. An attempt is made to distinguish these two cases:

1. If a value e has an imaginary part of very small magnitude, it's assumed that that is due solely to rounding error and Re(e) is returned;
2. If there's an imaginary part of larger magnitude, it's assumed that the value is truly non-real and undefined is returned.

The plot command doesn't care if some points are included that have undefined as one of their coordinates. It simply skips them. However, if a function has an undefined coordinate for all of its points, plot will issue a warning. To avoid these warnings, I remove those functions that always produce non-real values. There are 12 functions to start (as can be seen from the degree of W in x and y), and this filtering removes 4 of them.

So, here's the code that does all that: 

Digits:= 15: eps:= 1e-13: Ev:= evalhf: #All 3 may need to be adjusted.

#some testing values to decide real-or-complex issue; 99 is arbitrary:
T:= abc(b)/99 *~ {$1..99}:
 
(Px,Py):= (remove &< (w-> w~(T)={undefined}))~(
    (((w-> subs~&<(_P=~ w)) @ subs&<(V=~ z) @ [solve])~(W, V))( 
        z-> local e:= Ev(_P); `if`(abs(Im(e)) < eps, Re(e), undefined)
    )
)[]:
plot([`[]`~(Px, z->z, abc(rng))[], Py[]], abc(rng), scaling= constrained);

the

The above takes 0.516 seconds on my computer, and, because using evalhf makes the filtering and plotting nearly instantaneous, the vast majority of that time is spent in solve.

There's a chance (version-dependent I guess) that the solve will return RootOf expressions. In that case, add the explicit option. The polynomial is bicubic in both x and y, so it definitely can be explicitly solved.

There's a chance that some spurious horizontal lines will appear on the plot. These can be corrected by adjusting Digits, eps, or both. I haven't investigated the origins of these spurious lines. There's a remote chance, I guess, that you'll need to change evalhf to evalf to get higher precision, although I haven't yet seen a case where that's actually needed.