Probably your system can only be solved numerically:

restart;
odeSystem := {diff(y1(x), x) = -x*y2(x)-(1+x)*y3(x), diff(y2(x), x) = -x*y1(x)-(1+x)*y4(x), diff(y3(x), x) = -x*y1(x)-(1+x)*y4(x)-5*x*cos((1/2)*x^2), diff(y4(x), x) = -x*y2(x)-(1+x)*y3(x)+5*x*sin((1/2)*x^2), y1(0) = 5, y2(0) = 1, y3(0) = -1, y4(0) = 0}:
systemSol := dsolve(odeSystem, numeric):
plots:-odeplot(systemSol,[[x,y1(x)],[x,y2(x)],[x,y3(x)],[x,y4(x)]], x=0..1, color=[red,blue,green,gold], legend=[y1(x),y2(x),y3(x),y4(x)]); 

Here is another way to solve the problem using the custom procedure named  Circumcircle .  We can get both general formulas for the center and radius of the circumscribed circle of an arbitrary triangle, and specific results for a given triangle. We find the center of the circumscribed circle as the point of intersection of the two perpendicular bisectors to the sides of the triangle:

Download Circumcircle.mw

restart;
series(x/sin(x) + x, x = 0, 2);
convert(%, polynom);

                                                 1+x+O(x^2)
                                                       1+x

restart;
with(LinearAlgebra):
A := Matrix([[2, 3, 5], [1, 2, 7]]); 
P1:=ColumnSpace(A);
P2:=RowSpace(A);
 
B := Matrix([[6, 4, 2], [3, 2, 1]]); 
P3:=ColumnSpace(B);
P4:=RowSpace(B);

plots:-arrow(P1, width=0.03, color=red);
plot3d(t1*P2[1]+t2*P2[2], t1=-2..2, t2=-2..2);
plot([(convert(P3[],list)*~t)[],t=-3..3], scaling=constrained);
plots:-spacecurve(convert(P4[],list)*~t, t=-3..3, color=red, thickness=2);

Plotting.mw

The first of the conditions  f(x) = f(3 - x) means that the graph  y = f(x)  is symmetrical with respect to the straight line  x = 3/2 . The second condition simply sets the area under the graph  f(x)  on the segment  x = 1 .. 2 . But we do not know what values the function takes on the segment  x = 2 .. 3 . Therefore, it is impossible to calculate the last integral.
There may be a typo in the condition. If instead of the limits of integration  x = 2 .. 3  we take the segment  x = 1 .. 2, then the integral can be easily calculated using the symmetry and the second condition on the function  f(x) .

From wiki  https://en.wikipedia.org/wiki/Fermat_point :

1. Construct an equilateral triangle on each of two arbitrarily chosen sides of the given triangle.
2. Draw a line from each new vertex to the opposite vertex of the original triangle.
3. The two lines intersect at the Fermat point.

The following procedure finds the coordinates of the Fermat point in the triangle ABC. The formal parameters of the procedure are the vertices of the triangle, which are given as vectors and must be traversed counterclockwise.

restart;
FermatPoint:=proc(A::Vector,B::Vector,C::Vector)
local M, A1, B1, C1, Line;
M:=<0,1; -1,0>; # -90 degrees rotation matrix 
C1:=(A+B)/2+M.(B-A)*sqrt(3)/2;
A1:=(B+C)/2+M.(C-B)*sqrt(3)/2;
B1:=(C+A)/2+M.(A-C)*sqrt(3)/2;
Line:=(X,Y)->(y-X[2])*(Y[1]-X[1])-(x-X[1])*(Y[2]-X[2]); # Equation of a straight line through 2 points
solve({Line(A,A1),Line(B,B1)}, {x,y});
eval([x,y],%);
end proc:

Example:

FermatPoint(<0,0>,<7,0>,<3,5>);

Note: The method above applies if the largest angle in the triangle is less than 120 degrees. Otherwise, the Fermat point coincides with the vertex of the obtuse angle.

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