There are three dichotomies involved here, not just the two mentioned by the other respondents, and you are confusing them:

1. Standard vs. Classic GUI interface
2. 2D Input vs. 1D input (aka Maple Input)
3. Document mode vs. Worksheet mode

The red input that you've mentioned a few times is associated with 1D input. (You can change the color if you want.) You are confusing this with the other two dichotomies.

1. Standard vs. Classic: Classic is more robust and usually delivers the output faster. This has nothing to do with the speed of the computations performed by the kernel, only with the output---already computed---being formatted for your screen. There are a huge number of commands and command options (mostly related to plotting, graphics, tabular layouts) that will not work at all in Classic. Plus, you'll be restricted to a 32-bit kernel, which seems like a major setback.

2. 2D Input vs. 1D input (aka Maple Input): The 2D Input is to me nearly unreadable (due to a poorly spaced default font) and untypeable (due to constantly shifting context and the need to use arrow keys or the mouse). It is rife with bugs. Bugs are difficult to find because they may be due to character position rather than the character itself. Automatic conversion of 2D Input to 1D plaintext form results in a unreadable "brick". Despite what the other correspondents have said, there are a huge number of great syntax innovations that were introduced in Maple 2018 and 2019 that won't work at all with 2D Input. The error messages that you get for these will tend to make no sense; it's as if the 2D parser is totally unaware of that syntax. There are a moderate number of other forms of 1D syntax that work differently and erroneously when they are tried in 2D.

3. Document vs. Worksheet: I can't say much because I've never used Document.

So, I work exclusively with Standard GUI interface, 1D input (aka Maple Input), Worksheet mode.

OrbitPartition:= module()
option
	`Conceptual author: emendes`,
	`Maple code author: Carl Love <carl.j.love@gmail.com> 2020-May-13`
;
# Naming notes:
#    Names pertaining to the 1st index have "I" in the name: Is, nIs, PairsI, AllI, etc.
#    Names pertaining to the 2nd index have "J" in the name: Js, etc.
#    S always represents a tuple_size-combination of parms.

# Option cache or option remember?
#    No procedure with an 'S' argument in this module should be given option remember because
#    it takes too much memory and thus negates all the benefits of using Iterator. If such a 
#    procedure will be called again, give it option cache.
local
	#All module-level locals other than procedures
	#---------------------------------------------
	#module locals that are copied from parameters to ModuleApply:
	nIs::posint, #number of distinct 1st indices
	degree::nonnegint, #polynomial degree
	permutations_by_index::[table, table],
	CondTables::record,

	#module locals used in condition processing:
	Is::set, #set of first indices
	SubsetsI::set(set(posint)), #nonempty proper subsets of Is
	FullDegJ::posint, #critical value for condition FullDeg

 	#other module locals:
	parms::set([posint, nonnegint]), #fundamental set of pairs of indices
	nJs::posint, #number of distinct 2nd indices 
	Combos::Iterator,
	Setup::truefalse:= false, #Has problem been initialized?

	#Exclusion conditions:
	#---------------------
	#Build a table of sets of the 2nd index of S grouped by the 1st index:
	JsbyI:= proc(S)
	option cache, threadsafe;
		subs(
			_C= (op~)~(2, ListTools:-Classify(IJ-> IJ[1], S)),
			proc(i) option cache(nIs), threadsafe; `if`(assigned(_C[i]), _C[i], {}) end proc
		)
	end proc,
		
	#Make CondTables into a remember-table procedure:
     CT:= proc(K)
	option remember;
	     subs(_T= eval(CondTables[K]), proc(k) option remember; _T[k] end proc)
	end proc,
	
	Js:= proc(S) option cache; op~(2,S) end proc,

	Condition:= table([
		"AllI"=	 (S-> op~(1,S) = Is),
		"FullDeg"= (S-> max(Js(S)) >= FullDegJ),
	   	"ValidJ"=	 (S-> CT(':-ValidJ')~({Js(S)[]}) = Is),
	   	"FullDim"= (	 
	   		S-> not ormap(II-> `union`(JsbyI(S)~(II)[]) subset CT(':-FullDim')(II), SubsetsI)
   		)   		
   	]), 
   	#Conditions to check: The order is for efficiency, but that's just based on guesses.
	Conditions:= ["AllI", "FullDeg", "ValidJ", "FullDim"], 
			
	#Decide whether a permutation's orbit is allowed to be represented:
	AllowOrbit:= S-> andmap(k-> Condition[k](S), Conditions),	
	#===============================
	#End of exclusion-condition code

	#Set permutations to identity for unmentioned indices. Form Cartesian
	#product permutation of the permutations on the individual indices.
	T:= proc(T::table, ij) option remember; `if`(assigned(T[ij]), T[ij], ij) end proc,

	#This is the extension of permutation T to a permutation of tuple-size-combinations
	#of parms.
	Symm:= S-> map(IJ-> T~(permutations_by_index, IJ), S),

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

	#Process one chunk of combinations,
	#regardless of whether using sequential or multithreaded:
	DoTask:= proc(rn::[posint,posint])  #starting Iterator rank and number of iterates
	option threadsafe;
	local
		#See extended example on help page ?Iterator,RevolvingDoorCombination. 
		Combo_:= Object(Combos, 'rank'= rn[1]),
		has:= ModuleIterator(Combo_)[1], #We don't need the 'get' procedure; only need 'has'.
		C:= output(Combo_) #the Iteraror's Array output.
	;       
		(to rn[2] while has() do Orbit(parms[[seq(C+~1)]]) od)
	end proc,
		
	ModuleApply:= proc(
		tuple_size::posint, #tuple_size=1 => module setup only
		nIs::posint,
		degree::nonnegint,
		permutations_by_index::[table, table],
		CondTables::
			record(
				ValidJ::table,
				FullDim::And(table, set([set]) &under {entries})
			),
		#whether to use sequential or multithreaded (parallel) code:
		{sequential::truefalse:= false},
		#number of task "chunks", regardless of 'sequential':
		{numtasks::posint:= kernelopts('numcpus')},
		{nocompile::truefalse:= false} #Compile the Iterator?
	)
		#initializations:
		#----------------
		#Equate module locals to this procedure's parameters of the same name, as needed:
		Setup:= false;
		thismodule:-nIs:= nIs;
		thismodule:-degree:= degree;
		thismodule:-permutations_by_index:= permutations_by_index;
		thismodule:-CondTables:= CondTables;
		
		Init(2*min(numtasks, kernelopts(numcpus)));
		if tuple_size=1 then return fi; #problem setup only
		
		Combos:= Iterator:-RevolvingDoorCombination(
			nops(parms), tuple_size, ':-compile'= not nocompile
		);
		
		#Iterations:
		for local iter to 2 do #kludge to workaround error on 1st run after restart
			try
				return
					#Iterate over all "chunks" of combinations. (Chunks are
					#task-sized subsets of combinations.)
					`if`(sequential, map, Threads:-Map['tasksize'= 1])(
						DoTask,			
						{Iterator:-SplitRanks(Number(Combos), ':-numtasks'= numtasks)[]}
					)
			catch:
				if iter=1 then
					printf(
						"Error trapped, "
						"so efficiency measures are invalid for this run.\n"
					)
				else
					error
				fi
			end try
		od
	end proc,

	Init:= proc(cachesize::posint)
	local P, nterms:= (nIs,deg)-> binomial(nIs+deg, deg);
		Is:= {$1..nIs}; #1st indices
		FullDegJ:= nterms(nIs, degree-1);
		local Js:= {$0...(nJs:= nterms(nIs, degree))-1}, i, j;
		parms:= {seq(seq([i,j], i= Is), j= Js)};
		SubsetsI:= {seq(combinat:-choose(Is, i)[], i= 1..nIs-1)};
		forget~([for P in thismodule do P od])[];
		for P in [thismodule:-Js, JsbyI] do
			try 
				Cache:-Resize(cachesize, P)
			catch "procedure should have a cache":
				Cache(cachesize, procedure= P)
			end try
		od;		
		
		#Validate permutations:
		local ind, ent;
		for local ij,PT in permutations_by_index do
			if
				nops((ind:= {indices(PT, 'nolist')})) <>
				nops((ent:= {entries(PT, 'nolist')}))
				or not ind union ent subset op~(ij, parms)
			then
				printf("Index %d = %a.", ij, [Is,Js][ij]);
				error "%-1 table is not a permutation of its index", ij
			fi
		od;
		
		#Validate CondTables:
		if {indices(CondTables:-ValidJ, 'nolist')} <> Js then
			error "indices of ValidJ must be %1", Js
		fi;
		if not {entries(CondTables:-ValidJ, 'nolist')} subset Is then
			error "entries of ValidJ must be *sequences* from %1", Is
		fi;
		if {indices(CondTables:-FullDim, 'nolist')} <> SubsetsI then
			error "indices of FullDim must be %1", PairsI
		fi;
		if not `union`(entries(CondTables:-FullDim, 'nolist')) subset Js then
			error "entries of FullDim must be subsets of %1", Js
		fi;

		Setup:= true;
		return
	end proc,

	ModuleLoad:= proc() Setup:= false; return end proc
;
export
	CondCheck:= proc(S::set([posint, nonnegint]))
		if not Setup then error "first, call OrbitPartition to setup problem" fi;
		if not S::set([integer[1..nIs], integer[0..nJs-1]]) then
			error "index out range; nIs=%1, nJs=%2", nIs, nJs 
		fi;
		map(evalb@apply, Condition, S)
	end proc
;
	ModuleLoad()	
end module:

Download Symmetries.mw

A completely different solution is to extract the "y" values from a plot of the function. You don't need to look at the plot.

seq(op([1,1], plot(f, 0..1))[..,2]);

For finer control, this can be used with plot options numpoints, adaptive, and\or sample (see ?plot,options).

Plots can be combined with plots:-display. In the plot below, I used y1 > 0 for the first plot and y2 < 0 for the second so that they don't overlap.

restart:
initialset1 := {seq(seq([x1(0) = a, y1(0) = b], a = -2 .. 2), b = -2 .. 2)}:
DE1 := [diff(x1(t), t) = y1(t), diff(y1(t), t) = 1-x1(t)]:
p1:= DEtools:-phaseportrait(
    DE1, [x1, y1], t = -5 .. 5, initialset1, x1 = -2 .. 2, y1 = 0 .. 3, 
    stepsize = .1, color = green, numpoints = 600, thickness = 2, 
    linecolor = black
):
initialset2 := {seq(seq([x1(0) = a, y1(0) = b], a = -2 .. 2), b = -2 .. 2)}:
DE2 := [diff(x1(t), t) = -y1(t), diff(y1(t), t) = 1+x1(t)]:
p2:= DEtools:-phaseportrait(
    DE2, [x1, y1], t = -5 .. 5, initialset2, x1 = -2 .. 2, y1 = -3 .. 0, 
    stepsize = .1, color = green, numpoints = 600, thickness = 2, 
    linecolor = black
):
plots:-display(p1,p2);

There: That's substantial progress towards your goal. Now you just need to adjust the initial conditions so that the trajectories connect across the x-axis.

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.