You have some pretty complicated code. Here is a simple recursive procedure that returns all the fractions obtained in the nth step:

T:=proc(n)
local P;
option remember;
P:=proc(p)
local p1, p2;
p1:=numer(p); p2:=denom(p);
p1/(p1+p2), (p1+p2)/p2;
end proc;
if n=1 then return [1] else
[seq(P(T(n-1)[i]), i=1..nops(T(n-1)))] fi;
end proc:

Examples of use:
T(4);
sort(T(4));
sort([seq(T(i)[],i=1..4)]);  # Sorted list of all the fractions obtained in the first 4 steps 
                             

Addition. Similarly, we can write a procedure using the formula from wiki  q[i+1]=1/(2*floor(q[i])-q[i]+1) :

T1:=proc(n)
option remember;
if n=1 then return 1 else
1/(2*floor(T1(n-1))-T1(n-1)+1) fi;
end proc:

sort([seq(T1(n), n=1..100)]);

I rewrote the original code, simplifying it a bit. We can reduce the number of unknowns by 2 times if we write  C1 := <cos(s), sin(s))> , C2 := <cos(t), sin(t)> . Then the conditions  eq1  and  eq2  are satisfied automatically. We get that the vectors  C1  and  C2  either coincide or differ from each other by turning by an angle multiple of  90  degrees.
 

Edit.

Download Sol_edit.mw

I converted your 2D math input code to 1D code and replaced  `a__10 `  with a__10  and so on. Now everything works properly:
 

I advise you to immediately write code in 1D math input. It will be faster and much easier to track possible errors, as all characters are clearly visible.

Download coeff_new.mw

The syntax for this is  log[2](6) . Maple automatically translates this into the natural logarithm (with base e). For numerical calculation use the  evalf  command:

log[2](6);
evalf(%);

Here is another method that immediately provides expansions for  cos(6*t)  and  sin(6*t) . If you need to check it manually, then equate the results from the application of de Moivre's formula and Newton binomial:

evalc((cos(t)+I*sin(t))^6);
Re(%) assuming real;

Look at this question and the answer to it.

This is a well-known bug in some versions of Maple <= 2018. As a workaround, besides the depthlimit option, you need to increase  Digits :

restart;
kernelopts(version);
Digits:=20:
Optimization:-LPSolve(4*x1+5*x2+5*x3+2*x4,
                      {33*x1 + 49*x2 + 51*x3 + 22*x4 <= 120},
                      maximize, assume ={nonnegint}, depthlimit=10);

        Maple 2018.2, X86 64 WINDOWS, Oct 23 2018, Build ID 1356656             
                            [13, [x1 = 2, x2 = 1, x3 = 0, x4 = 0]]

It works:

Mat:=module()
    description "My Package";
    option package;
    uses Statistics;
    export LinModel; 
    LinModel :=proc(xliste::list, yliste::list,x::algebraic, $)
        return Fit(a*x+b,xliste,yliste,x)
        end proc; # LinModel 

    end module; #Mat

xliste:=[1,2,3]: yliste:=[2,4,5]:
Mat:-LinModel(xliste,yliste, x);
 

Only numerically:

f:=x->x^5 + 2*x^3:
g:=y->fsolve(y = x^5 + 2*x^3, x);  # This is the inverse function for f

Examples of use:

g(1);
plot([f, g], -3..3, -3..3, color=[red,blue], scaling=constrained);
evalf(Int(g, 0..2));
D(g)(2);
                                    

`/`(a, b); 
                       a/b

A simple procedure  LongSum  solves your problem. Formal parameters:  f  -  procedure for the common term, n0 - number of terms before the ellipsis,  n - position of the last term.

LongSum:=proc(f,n0,n)
`%+`(seq(f(k),k=1..n0))+`...`+f(n);
end proc:

Examples of use:
LongSum(n->x[n], 2, n);
LongSum(n->n/(n+1), 3, k);
LongSum(n->n/(n+1), 3, 100);

                           

If there is no any formula for the common term, then you can simply write (without a procedure) as in the example below:

a+b+c+`...`+f;

It is easy to modify this procedure to display not one, but several last terms or some intermediate terms, etc.

Download Check.mw

f2:=cos(y-Pi/3):
is(expand(f1/f2)=1);

                         true

Formally, the expressions  f1  and  f2  do not coincide for all values x and y. For example we get an undefined expression  0/0  for x=0, y=0 :
 

eval([numer(f1),denom(f1)],[x=0,y=0]);

                         [0, 0]
 

m := n->piecewise(type(n,even),0, 1):
f := proc(a,b)
if type(a*b,integer) then return m(a*b) else FAIL fi;
end:
 
f(a,b);
f(2,3);
f(1,3);
f(a,3);

sin(30*Pi/180);
30^`°`;