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

Try the  plottools:-getdata  command.