The problem is nearly impossible to accomplish without a graphical analysis, as below. Thus, I think the general problem is beyond the capability of computer algebra systems; although it may be possible to handle some simple cases programmatically.

Here is a graphical analysis of your problem:
 

Download SwitchIntegralOrder.mw

I hate this error message with a passion because it makes the naive user believe that they need to restart Maple, which is far more effort than is necessary, especially when they have multiple worksheets open. All that you need to do is

1. Save the worksheet;
2. Close the worksheet;
3. Reopen the worksheet.

Also, while it's theoretically possible that a firewall is the source of this problem, and that has actually happened in the past, the chance that that's actually true rather than being caused by a bug in Maple is close to 0.

Here is the OrbitPartition code parallelized with Threads and put into a module. This code is very conservative with memory, and it achieves a fairly impressive real-time speedup factor of n/2 where n/2 is the number of hyperthreaded cores with 2 threads per core. I tested on my two computers: one where n=4 and one where n=8.

Regarding my question on "as well as its counterpart(s)": I proved (by hand, not with Maple) that all 4 of your exclusion conditions (with the given current index permutation tables!) have this property: S satisfies the condition iff conds(S) does also.  In my terminology, violation of any condition necessarily excludes the entire orbit. This allowed me to simplify the code.

If a module only exports one procedure, say MyProc, then there's no need to export it. Instead, name the module MyProc and name the procedure ModuleApply. Then you call the procedure simply as MyProc(...). You'll see this in the code below.

OrbitPartition:= module()
option
	`Conceptual author: emendes`,
	`Maple code author: Carl Love <carl.j.love@gmail.com> 2020-May-7`
;
uses It= Iterator;
local 
	#Problem setup:
	#==============
	nIs:= 3, #number of distinct first indices 
	Is:= {$1..nIs}, #set of first indices

	#fundamental set of lists (indexed variables without a stem):
	i, j, parms:= {seq(seq([i,j], j= 0..9), i= Is)}, 

	permutations_by_index:= [
		table([2= 3, 3= 2]), #permutations of 1st index
		table([2= 3, 3= 2, 5= 6, 6= 5, 7= 9, 9= 7]) #... of 2nd index
	],

	#Exclusion conditions:
	#---------------------
	#Build a table of sets grouped by their 1st index:
	ClassifyI:= S-> ListTools:-Classify(x-> x[1], S),

	#Decide whether a permutation's orbit is allowed to be represented:
	AllowOrbit:= S->
		local C;
		Is = op~(1,S) and                                   #condition 1
		max(op~(2,S)) >= 4 and                              #condition 2
		nops({entries((op~)~(2,(C:= ClassifyI(S))))}) = nIs #condition 3
			and not                                        #condition 4:
		ormap(k-> op~(2,C[k]) subset [{0,1,4}, {0,2,7}, {0,3,9}][k], Is),
	#================
	#End of setup code

	#Set permutations to identity for unmentioned indices:
	T:= proc(T::table, k) option remember; `if`(assigned(T[k]), T[k], k) end proc,

	#Combine individual indices' permutations into a permutation function for lists:
	conds:= S-> map(x-> T~(permutations_by_index, x), S),

	#For any tuple of lists, compute its orbit under `conds`. Return a
	#representative of the orbit.
	Orbit:= proc(S)
	local r:= S, R;
		do R[r]:= () until assigned(R[(r:= conds(r))]);        
		{indices(R, 'nolist')}[1] #Return lexicographic min as representative.
	end proc,

	ModuleApply:= proc(tuple_size::posint, {sequential::truefalse:= false})
		local Combos:= It:-Combination(nops(parms), tuple_size);
		`if`(sequential, map, Threads:-Map['tasksize'= 1])(
			(rn::[posint,posint])->
				local 
					Combo_:= Object(Combos, 'rank'= rn[1]),
					has:= ModuleIterator(Combo_)[1], C:= output(Combo_), S
				;       
				(
					to rn[2] while has() do
						if AllowOrbit((S:= parms[[seq(C+~1)]])) then Orbit(S) fi
			 	 	od
				),
			{It:-SplitRanks(Number(Combos))[]}       
		)
	end proc
;
end module:

Download Symmetries.mw

You're very lucky that you're starting programming with Maple 2019 or later, because it has a very natural syntax for list building with "embedded loops", which makes building lists via tables usually unnecessary (although it's fine to use tables if you want). But don't acquire bad habits as you learn. The syntax L:= [op(L), x] is one of the worst habits.

See if you can understand this code. Feel free to ask questions about it:

#Returns remainder of a/b. Integer quotient is returned in q.
rm:= (a::nonnegint, b::And(posint, Not(1)), q::name)->
    a-b*(q:= trunc(a/b))
:
Reverse:= (L::list)-> #Reverses a list
   local k:= nops(L); L[[while k>0 do k-- od]]
:
#n is number to be converted; b is the base to convert to.
CNS:= (n::nonnegint, b::And(posint, Not(1)))->
    local r:= n; Reverse([do rm(r, b, 'r') until r=0])
:
CNS(12345, 10);
                               [1, 2, 3, 4, 5]

CNS(12345, 16);
                                [3, 0, 3, 9]

Of course, if this wasn't a basic programming exercise, I'd use the built-in irem instead of rm. And Reverse is equivalent to the stock command ListTools:-Reverse.

Maple makes all parts of this exercise extremely easy.

PythTripBinOp.mw

The following procedure does essentially the same thing Mmcdara's process, but it's much simpler. In particular, it doesn't require any manual intervention. This gives exactly the showstat format except that the statement numbers and the spaces that precede them are removed. So, if you prefer that format over "%P", or if you don't have Maple 2020 or later, then use this.

showstat_no_nums:= p->
    `debugger/printf`(
        "%s",
        StringTools:-RegSubs(
            "\n    "= "\n ", 
            StringTools:-RegSubs(
                "\n *[0-9]+"= "\n    ", debugopts('procpretty'= p)
            )
        )
    )
:
#example:
showstat_no_nums(showstat);
showstat := proc(p::{:-`::`, :-name, :-procedure, :-And(:-`module`,:-appliable)}, statnumoroverload::{:-`..`, :-integer}, statnum::{:-`..`, :-integer}, $)
local res;
    if _npassed = 0 then
        map(thisproc,stopat())
    else
        if _npassed = 1 then
            res := debugopts(('procpretty') = p)
        elif _npassed = 2 then
            res := debugopts(('procpretty') = [p, statnumoroverload])
        elif _npassed = 3 then
            res := debugopts(('procpretty') = [p, statnumoroverload, 
              statnum])
        end if;
        map2(`debugger/printf`,"\n%s",[res]);
        if not (procname::indexed and member('nonl',{op(procname)})) then
            `debugger/printf`("\n")
        end if
    end if;
    NULL
end proc

#Compare:
showstat(showstat);

showstat := proc(p::{:-`::`, :-name, :-procedure, :-And(:-`module`,:-appliable)}, statnumoroverload::{:-`..`, :-integer}, statnum::{:-`..`, :-integer}, $)
local res;
   1   if _npassed = 0 then
   2       map(thisproc,stopat())
       else
   3       if _npassed = 1 then
   4 >         res := debugopts(('procpretty') = p)
           elif _npassed = 2 then
   5           res := debugopts(('procpretty') = [p, statnumoroverload])
           elif _npassed = 3 then
   6           res := debugopts(('procpretty') = [p, statnumoroverload, 
                 statnum])
           end if;
   7       map2(`debugger/printf`,"\n%s",[res]);
   8       if not (procname::indexed and member('nonl',{op(procname)})) then
   9           `debugger/printf`("\n")
           end if
       end if;
  10   NULL
end proc

The error that you got is a very common problem. The solution is to solve the ODEs for their highest-order derivatives before passing to dsolve. If the solve returns multiple solutions, we must select one of them. In your case, the highest-order derivatives are x'' and y''. The solve command will (rightfully, but silently) object to there being no x'' or y'' in the first ODE. So, replace the left side of the first ODE with its derivative, which'll introduce the needed 2nd derivatives; since the right side is constant, just remove it (equivalent to setting it to 0). I renamed your pre-solved ODEs odes to emphasize that there will be multiple solutions. The solve command is

ode:= solve({odes}, diff~({x,y}(t), t$2))[1][]; #Use [1] or [2].

There are 2 solutions, both of which work in dsolve. I thought [1] was more interesting (after correcting your g to -9.81), but you'll need to decide which is physically appropriate. No other changes are needed to remove the error condition (but I'm not saying that your equations are correct!).

Regarding the physics: If we consider vertical to be the positive direction of y, as is usual, then g needs to be negative. Otherwise, I'll leave it to Rouben to address the physics.

plots:-odeplot(Sol, Loop, t= 0..3);

The most direct way to "produce a list with a for-loop" is

L:= [for i from 11 to 20 do i^2 od];
             L: = [121, 144, 169, 196, 225, 256, 289, 324, 361, 400]

The print command is intended for printing supplementary information. It should never be the way that you obtain the direct resuts of your code.

In the axis option, you can color the tick labels separately from the axes themselves, like this:

axis= [tickmarks= [3, subticks= 4, color= black], thickness= 0, color= gray]

Also, there are only 4 integers between 0 and 5, so I used subticks= 4. But if you really want to use 5, that's up to you.

You can make the gray as light as you want like this:

axis= [tickmarks= [3, subticks= 4, color= black], thickness= 0, color= COLOR(HSV, 0, 0, .85)]

The 3rd number after HSV is the fraction of white in the gray, so 0 is pure black and 1 is pure white. I usually like .85.

I generally set axesfont to [Helvetica, Bold, 8]. Give it a try; I find it easy it read without crowding the numbers.

I strongly suspect that you want to do the dot product of spherical vectors by (essentially) converting them to cartesian and then doing the ordinary dot product. Here is a short procedure for it:

#Deriving the Formula:
#---------------------
#Caution: The style of Maple's spherical coordinates wrt plotting is <rho, phi, theta>,
#whereas the style wrt VectorCalculus is <rho, theta, phi>, where in both cases
#phi (0..2*Pi) is the longitude and theta (0..Pi) is the latitude. The names of the 
#angular coordinates are not mathematically relevant and they may be switched to suit
#user preference, but their positions are relevant! The following conversion takes
#care of this switch.
#
factor(
    changecoords(      #This command uses the "plot" version of the coordinates, so
        changecoords(  #I relabel the vector positions.
            <<a1 | a3 | a2>>.<<b1, b3, b2>>, [a1,a3,a2], spherical
        ),
        [b1,b3,b2], spherical 
    )[1,1]
);
    a1 b1 (cos(a3) cos(b3) sin(a2) sin(b2)
       + sin(a3) sin(b3) sin(a2) sin(b2) + cos(a2) cos(b2))

#Implementing the Derived Formula:
#---------------------------------
`&dot_sph`:= (a::Vector(3), b::Vector(3))->
    #Assumptions: 1st entry of each vector is radius, 2nd is latitude (0..Pi), 
    #  and 3rd is longitude (0..2*Pi).
    simplify(
        a[1]*b[1]*(cos(a[3]-b[3])*sin(a[2])*sin(b[2]) + cos(a[2])*cos(b[2])),
        trig
    )  
: 
#Examples:
#---------
v1:= <3, Pi/10, Pi/4>:  v2:= <4, Pi/10, Pi/4>:
v1 &dot_sph v2;
                               12

<2, Pi/6, Pi/10> &dot_sph <3, 2*Pi/3, Pi/10>;  
                               0

<1, Pi/2, Pi/3> &dot_sph <2, Pi/4, Pi/6>;
                        1  (1/2)  (1/2)
                        - 3      2     
                        2              

I had two major revelations about what you're doing. (If I use any terms that you're unfamiliar with, please ask for definitions.)

1. The name stem (your alpha) is superfluous; all that matters are the names' indices, viewed as lists.
2. All that you're trying to do is partition the set of k-combinations of parms into the orbits of conds and select a representative from each orbit.

So, here's a procedure that uses that paradigm (points 1 & 2) to produce the same results as "something like"---but much faster, of course. This also corrects the issue that you mentioned regarding my prior update: The representative of each orbit is now its lexicographically minimal entry.

OrbitPartition:= proc(
    parms::And(set, listlist &under [op]),
    tuple_size::posint,
    permutations_by_index::list(table)
)
local 
    T:= proc(T,j) option remember; `if`(assigned(T[j]), T[j], j) end proc,
    conds:= varCoef-> map(x-> T~(permutations_by_index, x), varCoef),
    Done:= table(),
    Orbit:= proc(x)
    local r:= x, R:= table([x=()]);
        Done[x]:= ();
        do
            r:= conds(r);
            if assigned(R[r]) then break fi;
            R[r]:= (); Done[r]:= ()
        od;
        {indices(R, 'nolist')}[1]
    end proc,
    Orbits:= table(), C, i
;
    for C in Iterator:-Combination(nops(parms), tuple_size) do
        i:= parms[[seq(C+~1)]];
        if assigned(Done[i]) then next fi;
        Orbits[Orbit(i)]:= ()
    od;
    {indices(Orbits, 'nolist')}
end proc
:
parms:= {seq(seq([i,j], j= 0..9), i= 1..3)}:
T1:= table([2= 3, 3= 2, 5= 6, 6= 5, 7= 9, 9= 7]):
T2:= table([2= 3, 3= 2]):

newabc:= CodeTools:-Usage(OrbitPartition(parms, 7, [T2,T1])): 
memory used=8.16GiB, alloc change=103.77MiB, cpu time=108.64s, 
real time=63.40s, gc time=59.34s

nops(newabc);
                            1018628

If need be, the memory usage of the above can be reduced by eliminating Done, but this will increase the time by necessitating calling Orbit for every tuple.

To address your titular Question: Constructing and searching variable-sized sets (but not lists) is best handled by tables whose indices (aka keys) are the set elements and whose entries are the superfluous (). The code above has three such tables: Done, R, and Orbits.

There is a stock package/object MutableSet for doing the same thing. I haven't tested, but I doubt that it can beat the efficiency of my table-based method.  

Major syntax enhancements in Maple 2019 allow the same thing to be done without the tables and with greater efficiency. I'll put Maple 2019 code in a Reply.

The correct result can be obtained via a custom `evalf/...` procedure that doesn't need to consider the numerical-analytical intricacies of ln(1+x):

restart:
ln1:= x->  #ln(1+x)
   `if`(x::float and x > -1, evalf('procname'(x)), 'procname'(x)) 
:
`evalf/ln1`:= proc(X);
 local x:= evalf(X), oldDigits:= Digits;
    Digits:= max(Digits, Digits-op(2,x));
    local r:= evalf(ln(1+x));
    evalf[oldDigits](r)
end proc
:    
evalf[32](ln1(exp(-64)));
                                                -28
            1.6038108905486378529760870340137 10   

This works off the (correct) assumption that ln's own evalf procedure (see showstat(`evalf/ln`)) will work correctly provided that its argument (the 1+x) is evaluated with sufficient precision. Maple's whole network of evalf procedures is built this way.

Newton-Raphson is my preference. You can do it like this:

F := ((-288*(lambda + (2*k)/3 - 2/3)*sqrt(2)*lambda^2*arctan(sqrt(4*lambda/
    (1 - k) - 2)*sqrt(2)/2)/(1 - k)^2 - 144*(lambda + (2*k)/3 - 2/3)
    *Pi*lambda^2*sqrt(2)/(1 - k)^2 - 36*((lambda + (2*k)/3 - 2/3)*
    (4*lambda/(1 - k) - 2) + (10*lambda)/3 + (4*k)/3 - 4/3)*
    sqrt(4*lambda/(1 - k) - 2))*(1 - k)^2)/(256*lambda^2)
:
Digits:= 15:
N:= unapply(eval(lambda - F/diff(F, lambda), k= 1/10), lambda):
Newton:= proc(N, x0, tol::postive:= 10.^(1-Digits), maxiter::posint:= 99)
local r:= x0, old:= x0+1, k;
    for k to maxiter while abs((r-old)/r) >= tol do
        (r, old):= (evalf(N(r)), r);
        print(r);
    od;
    if k > maxiter then error "did not converge" fi;
    r
end proc
:    

Newton(N,1);
                       0.570668758044113
                       0.551790249418825
                       0.551586825212882
                       0.551586798435436
                       0.551586798435435
                       0.551586798435435

Of course, it's not necessary to print every iterate; I just thought that you might want that "to develop interest", as you say.

You'll want to use a table indexed by the tuples rather than an Array for tab (from "something like"). The entry pointed to by the tuple index is irrelevant; I usually use () for irrelevant entries. The reason for using a table is that the lookup time (determining whether something is an index) is essentially constant, regardless of the size of the table; whereas lookup in an Array is proportional to its size.

So this is my table-based variation on your "something like":

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

T1:= table([2= 3, 3= 2, 5= 6, 6= 5, 7= 9, 9= 7]):
T2:= table([2= 3, 3= 2]):
T:= proc(T::table, j) option remember; `if`(assigned(T[j]), T[j], j) end proc:

conds:= proc(varCoef::set(specindex(alpha)))
option remember; 
    map(x-> alpha[T~([T2,T1], [op(x)])[]], varCoef)
end proc:

for i in abc do
    if assigned(tab[conds(i)]) then 
        unassign('tab[conds(i)]');
        tab[i]:= ()
    fi
end do:

newabc:= {indices(tab, nolist)}:

This will give no significant time improvement for the 3-tuples and 4-tuples, but there's a factor-of-2 improvement for 5-tuples and a factor-of-4 improvement for 6-tuples. That factor will continue to increase exponentially for larger tuples because its competing against exponentially-lengthening Array lookups.

All six variables equal to 0 is a solution that satisfies all 16 equations. I'd guess that this is not the solution that you want.