Example:

A:=Matrix([[1, 0.00009], [2, 1]]):

map(fnormal,A, 10, 1e-4);    

       [ 1   0. ]

       [ 2   1  ]

Of course, the constant a should not be an integer, otherwise the denominator would be 0. Unfortunately, Maple does not convert the sum of trigonometric functions in the product, so the denominator should be simplified by hand:

simplify(algsubs(1-cos(2*Pi*a)=2*sin(Pi*a)^2,-2*Pi*sin(Pi*a)/(-1+cos(2*Pi*a))));

                                                             Pi/sin(Pi*a)

The parameter  a  affects only the size of the curve and not its form. Therefore, for the plotting you can take a = 1.

To find the parametric equation of the curve make the substitution y=x*(tan(t))^3, substitute in the equation of the curve. First, express  x  through  t, then  y. Finally, you obtain

x=a*(cos(t))^3,  y=a*(sin(t))^3,   t=0..2*Pi

EM:=unapply([E01$a,E02$b,E03$c], a,b,c):

EM(1,3,2);

EM(2,2,1);

evalf(Int(sqrt(x^3+1), x=1..2));

2.129861293

I think I understand what you need. Use a double loop as follows:

L:=[0,1,2]:

M:=[]:

for i to nops(L) do

for j to nops(L) do

M:=[op(M), [L[i], L[j]]]:

od: od:

M;

If your polynomial is a quadratic form, then use the Sylvester criterion. For an arbitrary polynomial  the appropriate command in Maple is unknown to me.

Repetitive elements can easily remove from the list, making the set from the list.

L:=[0, 0, 1, 1, 2, 2]:

M:={op(L)}:

combinat[choose](M,2);

{{0,1},{0,2},{1,2}}

a:=Matrix([[4,3,4],[7,7,3],[7,3,1]]):

b:=Matrix([[6,3,5],[9,8,0],[8,8,1]]):

L:=[a,b,b,a]:

add(LinearAlgebra[DiagonalMatrix]([0$(2*i-2),L[i],0$(8-2*i)]),i=1..nops(L));

r:=25:

Sol:=dsolve({diff(x(t),t)=-10*x(t)+10*y(t), diff(y(t),t)=r*x(t)-y(t)-x(t)*z(t), diff(z(t),t)=(-8/3)*z(t)+x(t)*y(t), x(0)=-1, y(0)=-1, z(0)=1}, numeric):

plots[odeplot](Sol,[x(t), y(t), z(t)], 0..30, axes=normal, numpoints=10000);

It turns out a very interesting picture:

See  plots[matrixplot]  command.

f := (x, y) -> ((2*x^2+y^2*(x^2*(-2+sqrt(1/x^4+4/(x^2*y^2)+4/y^4+4))-1))/(4*x^2-2*y^2+4))^(1/2) ;

g := unapply(diff(f(x, y), x), x,y);

plot([g(.1, y), g(1, y), g(3, y)], y = 0 .. 10);

A:=Matrix([[2,2,8,5], [6,3,4,9], [5,5,7,4], [2,1,3,2]]);

interface(rtablesize=100):

B:=LinearAlgebra[DiagonalMatrix]([A, -A, A]);

Instead of L in fsolve command after comma write L=0..infinity

You have a rather complicated transcendental equation. Solve command does not solve such equations. Replace solve command by fsolve command.