I don't fully understand your request, but surely a single remove command can do the job in its entirety without any need of the indices k of elements removed and without any need to consider the changing size of the list.

Why does conds return a list, not a set? Why does this code need lists at all?

Do you want to remove all triples from abc that contain any member of conds(...), i.e., for which there's nonempty intersection? Then do this:

parms:= {seq(seq(alpha[i,j], j= 0..9),i= 1..3)}: #set
abc:= combinat:-choose(parms,3):

T1:= table([2= 3, 3= 2, 5= 6, 6= 5, 7= 9, 9= 7]):
T2:= table([2= 3, 3= 2]):
T:= (T::table, j)-> `if`(assigned(T[j]), T[j], j):

conds:= (varCoef::set(specindex(alpha)))-> 
    map(x-> `?[]`(alpha, zip(T, [T2,T1], [op(x)])), varCoef)
:
cond:= conds(abc[1]):
remove(s-> s intersect cond <> {}, abc); 

The phenomenon is well explained by VV. I just want to point out that there's a fundamental contradiction in using Digits:= n for n > 15 with evalhf. Thus, it's no surprise that nonsense results. The command evalhf works at a fixed precision that's often call double precision in other languages. It follows very precise and standardized rules that are set down by IEEE and followed by nearly all processor manufacturers.

Not only is your code unnecessarily lengthy, it is also inefficient timewise and unadaptable to different graph orders.

To partition a list, set, etc., into equivalence classes, use ListTools:-Classify, like this:

GT:= GraphTheory: LT:= ListTools:
TG:= LT:-Classify(
    GT:-Diameter, 
    [GT:-NonIsomorphicGraphs](
        7, 
        output= graphs, 
        outputform= graph, 
        restrictto= connected
    )
):

At this point, TG[k] is the set of all graphs from the list of diameter k; in other words, it's the equivalence class. The equivalence classes can be counted by

TGn:= nops~(TG);
             TGn := table([1 = 1, 2 = 373, 3 = 387, 4 = 82, 5 = 9, 6 = 1])

The kind of fine type discrimination that you're describing can be done as follows: Let T1 and T2 be any individual types or sets of types. Then the type "all things that are T1 but not T2" is constructed as And(T1, Not(T2)). I can't imagine any situation such as you're describing that can't be handled by this mechanism.

The :: operator is very convenient for type checking. For example,

if x::And(algebraic, Not({name, constant, numeric}) then ....

The conditional expression above is equivalent to type(x, And(...)).

Here is your code with some small improvements that I'd recommend, including the correction to your original problem. Now, each intermediary variable is assigned before it is referenced. (That's usually a good way to go, although it's not mandatory in a "symbolic" language such as Maple.)

restart:
Digits:= 15: #becuase you have 15-digit constants in the code.
tmin:=10: #°C 	
tmax:=70: #°C 	
Tmin:= evalf[5](convert(tmin, temperature, Celsius, kelvin)):	
Tmax:= evalf[5](convert(tmax, temperature, Celsius, kelvin)):
alfa:= 0.3:		
i:= 0:
for T from Tmin by 2 to Tmax+1.25 do
    a12:= 2305.28444347652 - 9.14490843016421*T + 0.00680052257590234*T^2:
    a21:= -6665.24838284836 + 46.0897018087247*T - 0.0694991633494123*T^2:
    x2:= 1-x1:
    tau12:= a12/T:
    tau21:= a21/T:
    G12:= exp(-alfa*tau12): 
    G21:= exp(-alfa*tau21):
    lng1:= x2^2*(tau21*(G21/(x1+x2*G21))^2+tau12*(G12/((x2+x1*G12)^2))):
    lng2:= x1^2*(tau12*(G12/(x2+x1*G12))^2+tau21*(G21/((x1+x2*G21)^2))):
    lnga1:= subs(x1=xa1,lng1):
    lngb1:= subs(x1=xb1,lng1):
    lnga2:= subs(x1=xa1,lng2):
    lngb2:= subs(x1=xb1,lng2):
    r1:= lnga1+ln(xa1)=lngb1+ln(xb1):
    r2:= lnga2+ln(1-xa1)=lngb2+ln(1-xb1):
    r:= fsolve({r1,r2}, {xa1= 0..0.22, xb1= 0.22..1}):
    i:= i+1;
    #Guard against the possibility that fsolve doesn't finish:
    if r::set(name= numeric) then
        xA1[i]:= eval(xa1, r): xB1[i]:= eval(xb1, r)
    else
        xA1[i]:= 'unknown'; xB1[i]:= 'unknown'
    fi;
    print(i, T, xA1[i], xB1[i]);
od:

I present this Answer because I think that the while the Iterator package is definitely worthwhile using for this (it is very efficient both in time and memory), its syntax is quite unusual. Here is the code to generate all 6-combinations of a 9-set:

S:= {'rand'()$9}; #any set of size 9
#Generate all 6-combinations:
for C in Iterator:-Combination(nops(S), 6) do
     S[[seq(C+~1)]]
od;

Note the double square brackets. Note that C is a 1-D Array (with lower index 1) of 0-relative indices. Note that 1 is added to each element of C. Note that the iterator is passed no information about the base set, S, other than its size.

If you get an error related to the compiler, rerun with the option compile= false to the Combination command.

There are two main ways that Iterators can be run:

1. The index variable, which is always a 1-D Array with lower index 1, is the index variable in an "in-do" or "for-in-do" loop.
2. The index variable is the index variable of seq, add, or mul.

If I enter all the needed multiplication signs to your input (thus,

p:= -(729*beta*(1/2*(-1/9*(-5/27*beta-1/9)*lambda^6-1/9*(1/9*beta^2*TT+
(10/9*TE+2/3*TT)*beta-1/9*TE)*lambda^4+5/27*(2/5*TT*(3/2*TT+TE)*beta+TE*
(TE-2/5*TT))*beta*lambda^2-1/9*TE*beta^2*TT*(-TT+TE))*p1(m,t)^2+(beta*
TT-1/3*lambda^2)^2*(-1/3*lambda^2+TE)^3*((p2(m,t))/(lambda^2-3*TE)+3/2*(
(lambda^2+TE)*p1(m,t)^2)/((lambda^2-3*TE)^3))))/((3*beta*TT-lambda^2)^3*
(3*TE-lambda^2)^2);

) and then do 

normal(p);,

I get something equivalent to, and fairly close in form to, your desired output.

You can also get a nice form from

simplify(p);

If your goal is to prove that the two forms are equivalent, then enter the second form (thus,

s:= -3*beta*p2(m,t)/(lambda^2-3*beta*TT)+3*beta^2*(5*lambda^2-3*beta*TT)*
p1(m,t)^2/(2*(lambda^2-3*beta*TT)^3);

), and just do

simplify(p-s);
      0

Modulo the possibility of a bug (which seems extremely unlikely for simple rational functions such as these), simplification to 0 is a solid proof of equivalence (over the complex numbers).

Like this:

T:= table([2= 3, 3= 2, 5= 6, 6= 5, 7= 9, 9= 7]):
subsindets(
    varA, 
    `?[]`~(A, `[]`~(anything, {indices(T,'nolist')})), 
    a-> B[op(1,a), T[op(2,a)]]
);

As a single call:

map(op~, [1,2], ans1)

If you're using Maple 2019 or later, you can do this:

A:= Array([4]);
a:= A[-1]: to 3 do A,= ++a od;

The ,= is the new "append" operator, which is especially useful for efficiently appending to Arrays. 

The ++ is a new operator that does exactly what it does in C++: ++a adds 1 to a, stores the result in a, and returns the updated value.

The A[-1] gets the last element of A (which is also the first element in this case). This is not new syntax.

You can use seq like this:

A:= Array([4]):
n:= 3:
u:= upperbound(A): 
a:= A[-1]: A(u+1..u+n):= <seq(a+i, i= 1..n)>;

To me, that's nowhere near as clear as the new syntax.

All of the above assumes that you want to append to the end of the Array. If instead you want to insert new elements somewhere other than the end, that's still possible to do, but the commands required are quite different. Let me know if that's what you want. 
 

It seems that the piecewise expressions returned by PDF confuse the numeric integrator's ability to make an accuracy estimate. The simplify command applied to the integral (in Int or int ​​​​​​form) will select the relevant branch of the piecewise, and the numeric integration will finish.

Your equation has its variable, q, under the square root. As you may recall from high-school algebra, solving such an equation involves first converting it to an equivalent polynomial equation. Some of the polynomial's solutions may not work in the original equation because the square root indicates only the positive branch. Because your equation also has the parameters c, m, and s, there is no single expression for q that will work for all values of the parameters. That's why back-substituting and simplifying doesn't give you 0.

There is, however, an elaborate process by which you can classify these conditional solutions and know what are the conditions that the parameters must satisfy for each to be valid. It uses the parametric option to the solve command. I show it in the attached worksheet.

Download SolveParametric.mw

Here are the highlights of the worksheet:

restart:
p:= (1-u)*(u-m)-q*(u+c)+u*(1+m-2*u)-s:
u:= (1+m-q-sqrt((1+m-q)^2-4*(m+c*q)))/2:
p:= simplify(p);

The command evala(Norm(p)) finds the minimal polynomial equivalent to p. In other words, it removes the radical. 

pp:= evala(Norm(p));
         2  2      2            2                  2          
pp := 4 c  q  - c m  q + 2 c m q  + 6 c m q + 2 c q  - 4 c q s

      3      2      2                  2              2        
   - m  + 2 m  q + m  s - 3 m q s + 2 q  s - c q + 2 m  + 2 m q

                      2        
   - 2 m s - 3 q s + s  - m + s

I will assume that you're interested only in real solutions. Thus, the expression under the radical must be nonnegative. The parametric option to solve finds conditions on the parameters (c, m, and s in this case) that make certain solution branches correct. It only works on polynomials.

The following command takes about 30 seconds.

pwsol:= solve({pp, q^2 + (-4*c-2*m-2)*q + (m-1)^2 >= 0}, q, parametric, real);
nops(pwsol);
                              1731

That means that there are 866 cases! I list them in the worksheet.

My preferred alternative to Acer's 'SF' is :-SF. The :- prefix unmistakably marks it as global (to a Maple-savvy reader). Preben's global SF also unmistakably marks it as global, of course, but I prefer the :- because it's targetted to the specific instance of the variable and not other possible occurences in the same procedure.

':-SF' also works.

First, data and d need to be Vectors, not Matrixs!

(`+`@seq/numelems)~(
    map2(
        `?[]`, 
           d, 
        `[]`~(remove(`=`, [ListTools:-Split(j-> data[j]>0.01, [$1..numelems(data)])], []))
    )
); 
             [22/3, 51/2]

Here's one way:

{seq(seq(seq(seq((a*x0+b*x1)*(c*y0+d*y1), a= {0,1}), b= {0,1}), c= {0,1}), d= {0,1})};