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.

The only real solution of your system is a "singular" solution at x=0, f=57/40, which is essentially what fsolve is reporting. (You'll see what I mean by "singular" in the worksheet.) The plot you showed must involve different values of the parameters. The above solution is easy to get exactly, even by hand. With significantly more effort, I verify that this is the only real solution.
 

Download SingularSoln.mw

If you want to simplify an expression under the assumption that all of its variables are real, use evalc. For example,

evalc(temp1b)

will remove all the overbars, because conjugation is the identity function on real expressions.

That depends on whether you know, numerically, where the "top" is. If you know that the top is at 10, for example, you could do

plot([[5,0], [5,10]])

This is how the same thing can be achieved with the updating append operator ,=. The benefit of this operator is that you don't need to keep track of or increment the indices within the loop.

V:= Array(1..0): #empty
for k to 10 do
    V,= k/10.
od:
rtable_options(V, subtype= Vector[row]);
V;
[0.1000000000, 0.2000000000, 0.3000000000, 0.4000000000, 0.5000000000, 
 0.6000000000, 0.7000000000, 0.8000000000, 0.9000000000, 1.000000000]

For reasons that I don't understand, you must start with an Array rather than a Vector to do this. But the rtable_options command converts it to a Vector inplace, so that's a constant-time trivial operation (i.e., regardless of the number of elements). Anyway, it's very likely that an Array will work just as well as a Vector for your purpose.

Like all of the basic syntax features added in the last three major releases (Maple 2018 - 2020), the ,= operator is only usable in 1D input (aka Maple Input).

The following test plot shows that the lines are laid down in the order that they appear in the code:

plot([[[0,0],[1,1]],[[1,0],[0,1]]], thickness= 20);

However, when I tested the same thing a very long time ago (nearly 20 years), the lines were laid down in the reverse order that they appeared in the code. Since you're using a very old version of Maple, you should do this simple test.

Modifying my previous Answer (to which Acer linked) for non-constant theta(t) is easy. Here it is. I've highlighted in green the parts that you might want to vary. This code uses embedded assignments (highlighted in violet), which were introduced in Maple 2019 (maybe 2018). If you're using earlier Maple, I can easily modify this. Let me know.

plots:-display( #needed to get polar axes
    plottools:-transform((t,r)-> r*~[cos,sin](t))( #transform to polar
        DEtools:-phaseportrait(
            [
                diff(r(t),t) = r(t)^2*(r(t)+sin(theta(t))), 
                diff(theta(t),t) = r(t)^3*(r(t)+theta(t)+cos(theta(t)))
            ],
            [theta,r](t), 
            t= 0..0.5, #independent variable range for all trajectories 
            (i0:= map2(`~`[`=`], [theta,r](0), `[]`~(Pi/8*~[$0..15], .75))), #inits
            linecolor= [seq(COLOR(HSV, .85*i/ni, .85, .7), i= 0..(ni:= nops(i0)-1))],
            thickness= 1, #trajectory thickness (min. is 0, not 1)
            method= rkf45, stepsize= 0.0005,
            color= COLOR(HSV, 0, 0, .85), #light-grey arrows
            arrowsize= .3
        )
    ), 
    axiscoordinates= polar, size= [800$2], 
    axis[2]= [tickmarks= [subticks= 0]], axis[1]= [gridlines= false]
);

Notice that the trajectories that start at Pi/4 <= theta(0) <= Pi/2 are much longer than the others.