An n x n matrix with a column of zeros (such as your first example) cannot possibly have rank n, regardless of the algebraic field used. The rank of both of your example matrices is 2.

Change [x(0)=1, y(0)=1] to [[x(0)=1, y(0)=1]]. The reason is that DEplot allows you to specify multiple sets of initial conditions in a single plot, so the syntax requires them to be subgrouped.

Also, if you haven't loaded package DEtools with with(DEtools), then change DEplot to DEtools:-DEplot. I strongly prefer this package:-member syntax instead of with when one doesn't need the whole package. 

I get this plot

You almost had it. As you figured, the fundamental command for this is combinat:-randperm. So, you permute the right sides of the equations, and then match them with the left sides. The matching is most easily done with =~ (the elementwise equation-building operator):

([f,g,h,i] =~ combinat:-randperm([f||(1..4)]))(x);
           [                  x          2         1]
           [f(x) = x, g(x) = 3 , h(x) = x , i(x) = -]
           [                                       x]

It can be done like this:

alpha:= [15, -5/2, 15/11, -15/16, 5/7]:
p:= 3:  #starting position
beta:= alpha[[seq(p..1, -1), p+1..-1]];

I'm not sure if this'll work in your older Maple. If it doesn't, let me know; there'll certainly be something similar that works.

This is probably obvious: The elements of alpha can be anything; they needn't be numbers, real or otherwise. 

Use ':-p' < ':-c'. The purpose of the :- is to make it refer to the global variable rather than the parameter. The purpose of the single quotes is to make it ignore any values that may be assigned to the globals. Unevaluation techniques (such as single quotes) usually don't work with procedure parameters, which is why you need to switch to the globals.

Suppose that your DataFrame is named DF. To remove all rows for which any cell is undefined, do

DF[andseq(DF[C] <>~ undefined, C= with(DF))];

Great Question, by the way; it's not at all "simplistic". It took me about an hour of experimentation and help-page research to figure out how to do it.

Here's a significantly faster procedure:

OEIS_A219954:= proc(n::posint)
option remember;
local L:= ilog2(n), k:= n-1;
   `if`(n=1, 0, thisproc(k) + `if`(2^L=n, 3^L - n/2, 3^add(Bits:-Split(k))))
end proc
:

Edit: Changed the procedure from an arrow operator to a proc ... end proc to accomodate 2D Input (and I hope that I don't live to regret it) or older Maple.

First, I want to caution you to make functionList a proper Maple list by using square brackets [ ] instead of curly braces { }. The braces make it a set, and you can't control its order. (This is the only reason to make it a list rather than a set.)

The procedure below will take any of your functions and return side-by-side 3D plots of its real and imaginary parts.

ReIm_plot:= proc(f::algebraic, Zr::(name=range(complexcons)))
local z:= lhs(Zr), R:= rhs(Zr), x, y, RI:= [Re,Im];
    plots:-display(
        <
            plot3d~(
                RI(eval(f, z= x+I*y)), (x,y)=~ map~(RI, R)[], 
                labels=~ `[]`~(RI[](z), RI(evaln(:-f)(z))),
                title=~ typeset~(["Real", "Imaginary"], " part of ", f),
                _rest
            )
        >^%T
    )
end proc
:

functionList:= [2 + z, z^2 - 3*z, -z^3 + 4];

You can do one function via

ReIm_plot(functionList[3], z= -1-I..1+I));

Or you can do them all at once via
print~(ReIm_plot~(functionList, z= -1-I..1+I)):  #Note the colon terminator
 

>  functionList:= [2 + z, z^2 - 3*z, -z^3 + 4]: >  ReIm_plot:= proc(f::algebraic, Zr::(name=range(complexcons))) local z:= lhs(Zr), R:= rhs(Zr), x, y, RI:= [Re,Im];     plots:-display(         <             plot3d~(                 RI(eval(f, z= x+I*y)), (x,y)=~ map~(RI, R)[],                 labels=~ `[]`~(RI[](z), RI(evaln(:-f)(z))),                 title=~ typeset~(["Real", "Imaginary"], " part of ", f),                 _rest             )         >^%T     ) end proc : >  print~(ReIm_plot~(functionList, z= -1-I..1+I)):                         >

Download ReIm_plot.mw

I think that you want cylindrical coordinates. But showing the discontinuity (division by 0) elegantly is more difficult in 3d than 2d, regardless of your coordinate system.

plot3d(
    [r, theta, 1/(r^2*sin(theta)^2)], r= 0..2, theta= -Pi..Pi,
    coords= cylindrical, shading= zhue, view= [-2..2, -2..2, 0..99]
);

Your 3rd example shown in your Question does not have the problem that you're asking about. It doesn't change the abs.

To get just the number, use

eval(phi1, S2[2])

Edit: corrected argument order.

Here is a procedure for it:

IsRectangular:= proc(L::anything, max::posint:= infinity)
local d, Op:= op@{op}, Nops:= nops, Op0:= curry(op, 0), t;
    if not L::':-list' then return false, 0 fi;
    for d to max do until 
        (t:= (Op0:= Op@eval(Op0)~)(L)) <> ':-list' or 
        nops({(Nops:= Op@eval(Nops)~)}(L)) <> 1
    ;
    d=max or max=infinity and not ':-list' in {t}, d     
end proc
:

You still have exp^(-3*x[i]). You need to remove the ^ from that.

It can be done by repeating the first point in the plot command, like this:

Matrice:= [[-2,-3],[-1,2],[3,4],[1,-2]]:
plot([Matrice[], Matrice[1]]);
 

You need to use hastype instead of has, like this:

if pdext_simplified::`+` then
    (expr, other):= 0, pdext_simplified;
    for k from -1 to 1 do
        (yes, other):= 
            selectremove(hastype, other, u(anything$2, identical(n+k)));
        expr:= expr + ``(yes)
    od
fi:
expr;

This produces an expression with subgroupings enforced by ``(...), which should be adequate for display purposes. The ``(...) can be removed by expand(expr), but this will return it to its original term order. It is possible to sort, but it's a bit more complicated.