Usually  return  is used only in procedures. In order to find the positions of all columns that have duplicate elements, you can simply save these positions as values of some indexed variable, and at the end convert these values to a set to exclude repetitions.

Example:

restart;
Dati:=LinearAlgebra:-RandomMatrix(5, 10, generator=1..10);
k:=0:
for j from 1 to 10 do
for i from 1 to 5 do
h:=Dati[i,j];
for x from 1 to 5 do        
if Dati[x,j] = h then
if x <> i then    k:=k+1; L[k]:=j;
end if; end if;
end do: end do: end do:
L:=convert(L, set);

Addition. Of course, this problem can be more effectively solved, for example, so

seq(`if`(numelems(Dati[..,k])<>nops(convert(Dati[..,k], set)), k, NULL), k=1..10);
                                        1, 3, 4, 5, 6, 9, 10

Example:

interface(rtablesize=infinity):
A:=LinearAlgebra:-RandomMatrix(30, 10, generator = 1 .. 9);
seq(convert(A[..,k], list), k=1..10);

In fact, your polynomial   r^2 – r*s + s^2  is a positive definite quadratic form. This means that the inequality   r^2 – r*s + s^2>0  holds for any values of the variables  r  and  s  that are not simultaneously equal to 0. Sylvester's criterion (see https://en.wikipedia.org/wiki/Sylvester%27s_criterion) checks the quadratic form of any number of variables whether it to be sign-definite. I did not find the appropriate procedure in Maple and therefore wrote my own procedure with the name Sylvester:

Sylvester:=proc(Q::polynom, Var::list)
local Coeff, n, A;
uses LinearAlgebra;
Coeff:=proc(P, Var, t)  # Extracts the coefficient in front of the monomial  t of a multivariate polynomial  P  of variables  Var
local L,H,i,k:
L:=[coeffs(P, Var, 'h')]: H:=[h]: k:=0:
for i from 1 to nops(H) do
if H[i]=t then k:=L[i] fi:
od:
k;
end proc;
n:=nops(Var);
A:=Matrix(n, (i,j)->`if`(i=j, Coeff(Q,Var,Var[i]*Var[j]), Coeff(Q,Var,Var[i]*Var[j])/2));
if `and`(seq(Determinant(A[1..k,1..k])>0, k=1..n)) then return
`Positive definite` elif
 `and`(seq(Determinant(-A[1..k,1..k])>0, k=1..n)) then return
`Negative definite` elif
 `and`(seq(Determinant(A[1..k,1..k])>=0, k=1..n)) then return
`Positive semidefinite` elif
`and`(seq(Determinant(-A[1..k,1..k])>=0, k=1..n)) then return 
`Negative semidefinite` else
`Sign-indefinite` fi;
end proc:

Examples of use.

Your example:

Sylvester(t^2-t*s+s^2, [t,s]);
                                                       Positive definite

More complicated examples:

Sylvester(2*x1^2-8*x1*x2-4*x1*x3+9*x2^2+19*x3^2, [x1,x2,x3]);
                                                        Positive definite

Sylvester(x1*x2+x2*x3-x1^2-x2^2-x3^2,[x1,x2,x3]);
                                                        Negative definite

restart;
a:=1000:  
b:=100:    
x[1]:=a-b:  k[1]:=b:
for i from 2 to 30 do
k[i]:=2*k[i-1];  x[i]:=2*x[i-1]-k[i];
od:
z=add(x[i], i=1..30);

                                                     z = -2040109466700

Addition. Of course, this problem is better solved using  rsolve  command. This allows us to obtain explicit formulas for the sequences  x(n)  and  k(n) :

restart;
a:=1000:  
b:=100:
k:=unapply(rsolve({k(n)=2*k(n-1), k(1)=b}, k(n)), n);    
x:=unapply(rsolve({x(n)=2*x(n-1)-k(n), x(1)=a-b}, x(n)), n);
z:=unapply(sum(x(n), n=1..k), k);
z(30);
z(1000);

restart;
local D:     
d := s*x*(E*K*q-K*r+K*sigma[1]+r*x)*
     (1-x*(E*K*q-K*r+K*sigma[1]+r*x)
     /(K*sigma[2]*L))/(K*sigma[2])+sigma[1]*x
     -x*(E*K*q-K*r+K*sigma[1]+r*x)/K:
L := ldegree(d,x):
P:=d/x^L:    
x^L * (add(``(simplify(coeff(P,x^k),size))*x^k, k=[3,2,1])+``(simplify(tcoeff(P,x),size)));
([A,B,C,D]=~ListTools:-Reverse(expand~([coeffs(op(2,%), x)])))[ ];

The issue is that  gcd  command only works for two polynomials. Therefore, the simplest solution for your example is to apply it 2 times:

P:=[x^2-7*x+10, x^2, x^2+2*x+1];
gcd(P[1], gcd(P[2], P[3]));
                                                                 1

If you have many polynomials, then you can use a simple user procedure gcdp :

gcdp:=proc(P::list(polynom))
if nops(P)=2 then return gcd(op(P)) else
gcd(P[1], gcdp(P[2..-1]));
fi;
end proc:

Example of use for your example:

gcdp(P);
                                                      1

We obtain an equation of not third but fourth power of  x :

collect(d, x);
                     

restart;
f1 := x^2+y^2+z^2 = 1;
f2 := x+y+z=1;
with(plottools):
with(plots):
Circle:=intersectplot(f1,f2, x = -1 .. 1, y = -1 .. 1, z = -1 .. 1, color=red, thickness=3, numpoints=5000):
S1 := implicitplot3d(f1, x = -1 .. 1, y = -1 .. 1, z = -1 .. 1, style = patchnogrid, color = blue):
S2 := implicitplot3d(f2, x = -1 .. 1, y = -1 .. 1, z = -1 .. 1, style = patchnogrid, color = gold, transparency=0.5):
display(S1,S2,Circle, scaling = constrained, axes = boxed);
sol:=eliminate({f1,f2}, z);
Sol:=solve(sol[2,1], {x(t),y(t)});
assign(sol[1,1]):
assign(Sol):
Eq:=[x,y,normal(z)]; # Parametric equation of the circle
spacecurve(Eq, t=-500..500, color=red, scaling = constrained, axes = boxed, numpoints=20000);  # Visual check
 

Addition: another way.

Another way to derive parametric equations of this circle is to apply the rotation matrix. Let  d  be the distance from the origin to the plane  x+y+z-1=0  and  r=sqrt(1-d^2) .  Then this circle can be obtained from the standard circle  [x=r*cos(t), y=r*sin(t), z=d] (t=0..2*Pi)  by two turns: about the y-axis by the angle  arccos(1/sqrt(3))  and about the z-axis by the angle  Pi/4 :

restart;
local gamma:
beta:=arccos(1/sqrt(3)): gamma:=Pi/4:
A:=<cos(gamma),-sin(gamma),0; sin(gamma),cos(gamma),0; 0,0,1>.<cos(beta),0,sin(beta); 0,1,0; -sin(beta),0,cos(beta)>;
d:=eval(abs(x+y+z-1)/sqrt(3), [x=0,y=0,z=0]);
r:=sqrt(1-d^2);
eq:=<r*cos(t),r*sin(t),d>;
Eq:=simplify(convert(A.eq, list)); # Parametric equations of this circle
plots:-spacecurve(Eq, t=0..2*Pi, color=red, scaling = constrained, axes = normal);

plots:-polarplot(2*sqrt(cos(2*t)), t=0..2*Pi, color=red);

or (in cartesian axes)

plots:-polarplot(2*sqrt(cos(2*t)), axiscoordinates = cartesian, color = "Black", numpoints=5000, scaling=constrained);

or

plot(2*sqrt(cos(2*t)), t=0..2*Pi, coords=polar, color = "Black", numpoints=5000, scaling=constrained);
 

for the curve  x=cos(3*t), y=sin(2*t), t = 0 .. 2*Pi . 

restart;
plot([cos(3*t), sin(2*t), t = 0 .. 2*Pi]);
Sol:=solve({cos(3*t1)=cos(3*t2), sin(2*t1)=sin(2*t2)}, explicit);
simplify([Sol]);
map(t->[cos(3*rhs(t[1])),sin(2*rhs(t[1]))], convert(%,set));
% minus {%[-1]};  # The final result

plots:-inequal (for nonlinear curves)  and  plots:- shadebetween  commands appear only in the latest versions of Maple. For older versions,  filledregions=true  option can be used (in new versions it works also):

f:= x-> -x^2+3*x+2:
l:= x-> 4*x:
m:= x-> 2*x-4:
plots:-display(plots:-implicitplot(max(y-f(x), y-l(x), m(x)-y, -y)<0,  x= -1.5..4.5, y= -1.5..5.5, coloring = ["LightCyan", white], filledregions = true, gridrefine=3), plot([f(x),l(x),m(x)], x= -1.5..4.5, y= -1.5..5.5, color=[red,green,blue], thickness=2), scaling=constrained);

The same scale for each axis has been selected (scaling=constrained option)
 

The procedure  Check  solves the problem for a square matrix.

Check:=proc(A::Matrix)
uses LinearAlgebra;
is(`*`(convert(Eigenvalues(A), list)[ ])=Determinant(A));
end proc: 

Example of use:

A:=LinearAlgebra:-RandomMatrix(3, generator=0..9);
Check(A);

If you have several matrices, then just apply this procedure to each of them.

evala(convert(c, radical));

                                                     0

1.  RootOf  is a special way of presenting a solution, which is not always useful.
2. You have a rather complicated non-linear equation with several parameters. I think that it can be solved only numerically for specific values of the parameters and with the specified initial condition.

Example of numerical solution:

restart;
Eq:=diff(y(x), x) = -A*y(x)/x+x^alpha*(C+(y(x)^alpha/x^alpha-C)^(-`gamma`))/y(x)^alpha:
Sol:=dsolve({eval(Eq, [A=1,C=2,alpha=0.5,`gamma`=3]), y(1)=2}, numeric);
eval(y(x), Sol(2));  # Calculation of the value y(2)
plots:-odeplot(Sol,[x,y(x)], x=1..5, view=[0..5,0..4]);  # The plot of y(x)

This gives a higher quality:

restart;
expr:=x^2+y^2-2*y+1:
plot3d([r*cos(t), r*sin(t), eval(expr, [x=r*cos(t), y=r*sin(t)])], r=0..2, t=0..2*Pi, axes=normal, scaling=constrained);

Addition. One more way.

We can regard this surface as the image of the disk  x^2+y^2<=4  under the mapping  (x,y) -> (x, y, x^2+y^2-2*y+1). This method is useful in some cases:

Disk:=plot3d([r*cos(t),r*sin(t),0],r=0..2,t=0..2*Pi, style=surface, color=khaki):
f:=plottools:-transform((x,y)->[x,y,x^2+y^2-2*y+1]):
plots:-display(Disk, f(Disk), axes=normal, scaling=constrained);  # The surface and its projection onto the xy-plane