testeq  command is not suitable for checking all the roots, you can use  eval  command.

Solution and check:

eqns := {x^3-4*x = y, y^3-4*y = x}:

vars := {x, y}:

solns := solve(eqns, vars, explicit):

solns1 := map(simplify, [solns]);  # Simplified solutions

seq(simplify(eval([lhs(eqns[1])-rhs(eqns[1]), lhs(eqns[2])-rhs(eqns[2])], solns1[i])), i = 1 .. 9); #Checking of all the solutions

You got the solutions in implicit form. Here is the complete list of the all solutions:

eqns := {x^3-4*x = y, y^3-4*y = x};
vars := {x, y};
solns := solve(eqns, vars);

map(allvalues, [solns]);

Example:

L := [1,2,3,4,3,4,4];

{ListTools[FindRepetitions](L)[]};  # repetitive elements

convert(L, set) minus %;  # unique elements

                                   L := [1, 2, 3, 4, 3, 4, 4]

                                                {3, 4}

                                                {1, 2}

Addition: With regard to the indices see  ListTools[Search]  and  ListTools[SearchAll]  commands

In menu  view  classic Maple 12 you can read   (and this works)  Zoom Factor-> 50%  ctrl+0,  75%  ctrl+1, and so on.

In standard this does not work.

I do not see any problems.

Example:

f:=x->piecewise(x<=-1,x^2-2, x>-1 and x<1,x, x>=1,-x^2+3);

plot(f, -2..2, thickness=2, discont, scaling=constrained);

P  procedure solves your problem.

restart; 

P:=proc(m::posint, n::posint) 

local roll, i, j, k, L, A;

roll:=rand(100..999);

for i to m do

for j to n do

k:=roll(); L:=convert(k,base,10);

if L[1]<>L[3] then A[i,j]:=`   ` else A[i,j]:=k fi;

od;  od;

Matrix(m,n, (i,j)->A[i,j]);

end proc:

Example of use:

interface(rtablesize=infinity):

P(14,20);

In standard Maple 12:

restart;

with(VectorCalculus):

SetCoordinates(cartesian[x, y, z]);

alias(u = u(t, x, y, z), v = v(t, x, y, z), w = w(t, x, y, z)); alias(eta = eta(t, x, y, z));

U := VectorField(`<,>`(u, v, w));

Divergence(U);

Jacobian(U);

map(Diff, U, t);

convert(%, Vector[row]);

Probably you are using an older version of the package,  which in fact does not have  DirectSearch:-SolveEquations  command. Here's a link to the new version  http://www.maplesoft.com/applications/view.aspx?SID=101333

The easiest way to solve your problem by using a special procedure. The procedure  RootsToTrig  is based on the formulas in your link. The construction  ``(...)  used to prevent premature cosine computation in some cases.

RootsToTrig:=proc(P::polynom)

local P1, A, B, C, Q, R, Discr, F, theta, x;

x:=op(indets(P));

P1:=P/coeff(P,x^3);

A:=coeff(P1,x^2); B:=coeff(P1,x); C:=coeff(P1,x,0);

Q:=simplify(3*B-A^2)/9; R:=simplify(9*A*B-27*C-2*A^3)/54; Discr:=Q^3+R^2;

if is(Discr>=0) then error "Should be discriminant<0" fi;

F:=simplify(R/sqrt(-Q^3));

theta:=arccos(F);

x[1]=2*sqrt(-Q)*cos(``(theta/3))-A/3, x[2]=2*sqrt(-Q)*cos(``(theta/3+2*Pi/3))-A/3, x[3]=2*sqrt(-Q)*cos(``(theta/3+4*Pi/3))-A/3;

end proc: 

Examples of uses (all angles are expressed in radians):

RootsToTrig(5*x^3-4*x+1);

RootsToTrig(5*x^3-4*x+2);

 Here is the version of the same procedure with degrees:

restart;

RootsToTrig1:=proc(P::polynom)

local P1, A, B, C, Q, R, Discr, F, theta, x;

x:=op(indets(P));

P1:=P/coeff(P,x^3);

A:=coeff(P1,x^2); B:=coeff(P1,x); C:=coeff(P1,x,0);

Q:=simplify(3*B-A^2)/9; R:=simplify(9*A*B-27*C-2*A^3)/54; Discr:=Q^3+R^2;

if is(Discr>=0) then error "Should be discriminant<0" fi;

F:=simplify(R/sqrt(-Q^3));

theta:=arccos(F);

x[1]=2*sqrt(-Q)*cos(theta/3*180*degrees/Pi)-A/3, x[2]=2*sqrt(-Q)*cos((theta/3+2*Pi/3)*180*degrees/Pi)-A/3, x[3]=2*sqrt(-Q)*cos((theta/3+4*Pi/3)*180*degrees/Pi)-A/3;

end proc: 

Example:

RootsToTrig1((x-1)*(x-2)*(x-3));

There is no  evaluate  command in Maple (there is  eval  command). Removal (or expansion) of parentheses in different meanings expand command carries out.

Examples:

restart;

L:=[x*(x+y), (x+y)^5, sin(3*x), cos(x+y), (x+y)/z, 2^(x-y), Int(f(x)+g(x),x)]:

for l in L do

expand(l);

od;

restart;

solve({x+y-3, x^2+y^2-5}, {x, y}):

a,b:= eval([x, y], %[1])[]:

a, b;

                    2, 1

Array( [seq(plot(p, x = -Pi .. Pi, -1.4 .. 1.4, scaling = constrained), p = [sin(x), cos(x)])]);
plots[display](%);

You can just copy and paste the sign (below I put this sign in this way)

         ▪

of the Word. In the Word menu  Insert->Symbols->Other Symbols . The number of the sign in Unicode is  25AA .

Addition:  long to not look for this sign, you can just type in the Word  25AA , select it and then  alt+x

This  is equivalent to the union of two plots  r=-sqrt(5-4*cos(theta))  and  r=sqrt(5-4*cos(theta)) :

plots[polarplot]([-sqrt(5-4*cos(theta)), sqrt(5-4*cos(theta))], theta = 0 .. 2*Pi, color = red, thickness = 2);

In this example, it is more convenient to use  surd  command.  surd(a,3)  means  the cubic root in the real domain. The results are displayed in a more compact form.

Compare:

restart;

de := diff(y(x), x) = 3*abs(y(x))^(2/3);

dsolve(de);

and

restart;

de := diff(y(x), x) = 3*surd(y(x)^2, 3);

dsolve(de);

The latter result can be more simplified to  x - y(x)^(1/3) +_C1 = 0 ,  but Maple does not want to do this.