Here is another short way:

x:= <1, 2, 3, 4>:
k:= <seq((x+~x[j])^%T,j=1..4)>;

Here is a short solution that uses the parametric equations of two lines, which OP wrote down in his comment. The position of the point   A  on the first line depends on the parameter  t , and the position of the point  B  on the second line depends on the parameter  s . We write out vectors  MA  and  MB and compose the system  with 2 unknowns  t  and  s  from the condition of collinarity of these vectors:

restart;
A:=<1+2*t,t,-3-2*t>:
B:=<1-s,s,3+2*s>:
M:=<1,-2,3>:
MA:=A-M: MB:=B-M:
Sol:=solve(convert(LinearAlgebra:-CrossProduct(MA, MB),list)=~[0,0,0]);
V:=eval(A-B,Sol)*11/9; # Direction vector of the desired line 
(<x,y,z>)=~M+t*V;  # Parametric equations of the desired line

We get 2 solutions for this problem:

restart;
P1:=[1,0,3]: P2:=[1,0,3]+[-1,1,2]:
Plane:=A*x+B*y+C*z+E:
M:=[1,0,2]:
Sys:={eval(Plane,[x,y,z]=~P1)=0,eval(Plane,[x,y,z]=~P2)=0,abs(eval(Plane,[x,y,z]=~M))/sqrt(A^2+B^2+C^2)=1/4, E>0};
Sol:=[solve(Sys)];
sort(collect(simplify(eval(Plane,Sol[1][1..3])/E*19),[x,y,z]))=0; # The plane 1
sort(collect(simplify(eval(Plane,Sol[2][1..3])/E*19),{x,y,z}))=0; # The plane 2

The point  M  itself lies on this straight line, so if the plane passes through this straight line, then the distance from  M  to the plane is  0 .

A simple procedure  Mean  finds the arithmetic mean of any number of numbers:

Mean:=proc() `+`(args)/nargs; end proc:

Example of use:

Mean(3.2, 5.4, 4.6, 7);

                                              5.050000000

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