When you solve a system of equations, parameters for different curves must have different names:

restart:
l1:=[57/5 + 20*t, -21/5 - 40*t, 20*t]:
l2:=[249/100 - 120*s, 471/100 - 80*s, 200*s]:
Sol:=solve(l1=~l2);
P:=eval(l1, Sol);
eval(l2, Sol); # Check

The issue is that the  coeffs  command does not return the coefficients of the polynomial in the order in which the polynomial was written. To set this order you can use the option  t  (see the code below). The simple user-defined procedure  coef(P, m, var)  returns the coefficient of a polynomial  P  in one or more variables  var  before the specified monomial  m :

restart:
Digits := 20:
unprotect(D);
G := 0.04361098108*x^2 + 0.4810001561*x*y + 1.326278064*y^2 - 0.7320831383*x - 2.656083763*y + 1 = 0;
f := (x, y) -> lhs(G);
coeffs(f(x,y),[x,y],'t'); # The order of the coefficients is broken
t;

coef:=proc(P,m,var)
local L1,L2,T,t;
L1:=[coeffs(P,var,'t')];
L2:=[t];
T:=table(L2=~L1);
`if`(m in L2,T[m],0);
end proc:

A,B,C,D,E,F:=seq(coef(f(x,y),t,[x,y]), t=[x^2,x*y,y^2,x,y,1]); The order of the coefficients is ok

If I understand your question correctly, then you have a directed segment defined in the  geometry  package (or in the  geom3d  package).
An example:

restart;
with(geometry):
point(A, 1, 1), point(B, 3, 2):
dsegment(AB, [A, B]):
L:=map(coordinates,DefinedAs(AB)); # List of segment ends coordinates
convert(L[2] - L[1], Vector);

You can use an element-wise solve:

restart;
solve~({{x + y = 2, y = 1, x > 0},{x^2 = 4, x < 0, y = 0}});
# Or
solve~([{x + y = 2, y = 1, x > 0},{x^2 = 4, x < 0, y = 0}]);

                                           {{x = -2, y = 0}, {x = 1, y = 1}}
                                           [{x = 1, y = 1}, {x = -2, y = 0}]

restart;
with(VectorCalculus):
SetCoordinates(cartesian[x, y, z]):
a := <5, 0, 0>;
b := <4, 3, 0>;
Norm(a);
Norm(b);
v := CrossProduct(a, b);
Norm(v);

Riemann Stieltjes integral can be calculated if it can be reduced to the usual Riemann integral. This can always be done for piecewise continuously differentiable functions. The integral itself must be written in Maple syntax (I use 1Dmath input):

restart;
alpha := x->floor(x):
f:=x->x^2+1:
int(f(x)*diff(alpha(x),x), x=0..3);

I don't know the reason for this bug. As a workaround, you can use the  subsindets  command:

subsindets(<sin(u),sin(u)>, sin(anything),t->2*sin(op(1,t)/2)*cos(op(1,t)/2));

From wiki: "In geometry, a set of points are said to be concyclic (or cocyclic) if they lie on a common circle". 
I arbitrarily took 4 points  [cos(-Pi/6),sin(-Pi/6)], [cos(Pi/4),sin(Pi/4)], [cos(3*Pi/4),sin(3*Pi/4)], [cos(-Pi/2),sin(-Pi/2)]  on the unit circle. A total of 4 solutions were found, in two of which the parabolas  F  and  G  coincide. In the figure - the first solution from the list  L1:

restart;
F:=a1*x^2+b1*x*y+c1*y^2+d1*x+e1*y+1=0:
G:=a2*x^2+b2*x*y+c2*y^2+d2*x+e2*y+1=0:
F1:=b1^2-4*a1*c1=0: G1:=b2^2-4*a2*c2=0:
P:=[[cos(-Pi/6),sin(-Pi/6)],[cos(Pi/4),sin(Pi/4)],[cos(3*Pi/4),sin(3*Pi/4)],[cos(-Pi/2),sin(-Pi/2)]]:
Sys:={seq(eval(F,[x,y]=~p),p=P),F1,seq(eval(G,[x,y]=~p),p=P),G1}:
L:=[solve(Sys, explicit)]:
L1:=evalf(L):
F:=eval(F,L1[1]);
G:=eval(G,L1[1]);
Parabols:=plots:-implicitplot([F,G], x=-3/2..2, y=-5/2..3/2, color=[red,blue], gridrefine=3):
Points:=plots:-pointplot(P, color="Red", symbol=solidcircle, symbolsize=14):
Circle:=plottools:-circle(color=cyan):
plots:-display(Parabols,Points,Circle, scaling=constrained, size=[500,500]);

                 
Addition. In fact, we have the only solution. In solutions  L[1]  and  L[4] , parabolas  F  and  G  differ only in the order of succession.

2_parabolas.mw     

restart;
expr:=exp(3*I*x - x):
evalc(expr);  
simplify(evalc(expr));
# Or 
factor(evalc(expr)); 

restart;
M:=<<seq(j+1,j=1..5)>|<seq(2*j, j=1..5)>>:
plot(M, style=pointline, color=red, symbol=solidcircle, symbolsize=12, view=[0..6,0..10], scaling=constrained);

It doesn't need loops. This is done directly from two vectors using the  seq  command:

restart;
<<x,seq(j+1,j=1..5)>|<y,seq(2*j, j=1..5)>>;

It can be made prettier using the  DocumentTools:-Tabulate  command:

M:=<<x,seq(j+1,j=1..5)>|<y,seq(2*j, j=1..5)>>:
DocumentTools:-Tabulate(M, width=10):

One intermediate step:

C1 := sqrt(r+1)*sqrt(1/(r-1)); C2 := -sqrt(r^2-1)/(r-1);
C11:=combine(C1) assuming r+1>0;
is(C11=C2) assuming r<1, r>0;

Your matrix  F  is only defined for values on the main diagonal and below. You can define its values above the main one as you like. For example, the code below defines a skew-symmetric matrix:

restart;
for i to 6 do
for j to i do
F[i, j] :=  g(i+1)/g(j);
 end do:
end do:
F:=Matrix(6,6,(i,j)->`if`(i>=j,F[i,j],-F[j,i]));

If I understand your question correctly, here are two simple ways:

restart;
interface(rtablesize=infinity):
Matrix([seq([i,ithprime(i)], i=1..99)]);
# Or
Matrix(99,2, (i,j)->`if`(j=1,i, ithprime(i)));

For the calculation, we used the evenness of the function  f(k)=(-1)^k/k * sin(k*x) , the value at zero  k=0  calculated  through the limit, and also used the assumption (assuming real)  for the parameter  x :

restart;
Sum((-1)^k/k * sin(k*x), k=-infinity..infinity)=limit((-1)^k/k * sin(k*x),k=0)+2*evalc(sum((-1)^k/k * sin(k*x), k=1..infinity)) assuming real;

The result obtained can be substantially simplified to a piecewise constant function (R11 in the code below). Unfortunately, Maple does not make such a complete simplification (does in the range  x=-Pi..Pi  only)  and the piecewise-function  R11  is found manually:

R:=x-2*arctan(sin(x)/(cos(x)+1));
R1:=normal(applyop(convert,[2,2,1],R, tan));
simplify(R1) assuming x>-Pi,x<Pi;
R11:=piecewise(x=(2*n-1)*Pi,(2*n-1)*Pi,x>-Pi+2*Pi*n and x<Pi+2*Pi*n,2*Pi*n); # assuming n::integer
plots:-display(plot(R, x=-5*Pi+0.01..5*Pi, color=red, thickness=2, discont),plot([seq([(2*n-1)*Pi,(2*n-1)*Pi],n=-2..2)],style=point, symbol=solidcircle, color=red), scaling=constrained);