To do something with the resulting closed curve, for example, calculate a double integral over the corresponding area or something else, you can first get an array from the coordinates of the points lying on this curve (below is the P1 array). We then actually approximate this curve with a polygon with vertices at those points. In the example, we find the perimeter of this area:

To calculate a double integral over this area of a function of 2 variables, you can use Green's formula.
 

Download contact-Problem_new1.mw

The problem has 2 solutions. It is clear that the coefficients  A, B, C, D  can only be found up to a constant non-zero factor. Therefore, at the end, we simply reduce the equations by this factor:

restart:
local D:
np:=<A,B,C>: nq:=<1,-1,0>: P:=A*x+B*y+C*z+D=0: M:=[1,0,0]: N:=[0,0,-1]:
Sys1:=eval(P,[x,y,z]=~M);
Sys2:=eval(P,[x,y,z]=~N);
Sys3:=(np.nq)/sqrt(np.np)/sqrt(nq.nq)=cos(Pi/4) assuming real;
Sol:=[solve({Sys1,Sys2,Sys3})];
eval(P,Sol[1])/B; # The first solution
simplify(eval(P,Sol[2])/B); # The second solution

restart:
plots:-implicitplot(x^2+(y-surd(x^2,3))^2 = 1, x=-3..3, y=-3..3, color=red, thickness=3, gridrefine=3);

                  ,     

I replaced  x^(2/3)  with  surd(x^2,3), because x^(2/3) for  x<0  returns complex values (see help on the  surd  command).

restart;
assume(d::'real', d>-infinity, d<infinity);
coulditbe(d=Pi);
coulditbe(d=infinity);
coulditbe(d=1+I);

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):