This is quite simply:

M1:=[-2,3,-4]:  M2:=[4,-1,3]:

A:=[-t+3,2*t-4,t+2]:  B:=[s+4,2*s,2*s-1]:

Dist1:=sqrt((M1[1]-A[1])^2+(M1[2]-A[2])^2+(M1[3]-A[3])^2)+sqrt((B[1]-A[1])^2+(B[2]-A[2])^2+(B[3]-A[3])^2)+sqrt((B[1]-M2[1])^2+(B[2]-M2[2])^2+(B[3]-M2[3])^2):

Dist2:=sqrt((M1[1]-B[1])^2+(M1[2]-B[2])^2+(M1[3]-B[3])^2)+sqrt((B[1]-A[1])^2+(B[2]-A[2])^2)+(B[3]-A[3])^2+sqrt((A[1]-M2[1])^2+(A[2]-M2[2])^2+(A[3]-M2[3])^2):

evalf(minimize(Dist1));  evalf(minimize(Dist2));

min(%, %%); 

                                                                         21.59170144

                                                                         16.24881940

                                                                         16.24881940

1) If you consider the function y=psi(x)  in the real range, then you should write surd(1-x, 3)  instead  (1-x)^(1/3). Else for  (1-x)^(1/3)  if  x>1 , you get complex values.

2) By the code  q:=y->fsolve(y=psi(x), x)  you are looking the function inverse to the function y=psi(x) . For correct solution of the problem the function y=psi(x) should be monotonic in the corresponding range. Such ranges simpliest can be found using the plot of the function. Function  y=psi(x)  in the entire domain does not have an inverse function.

V:=1.8e-5:  R:=8.314:  T:=298:  k:=0.841:  G:=1e6:

psi:=unapply((ln(x)-x+k*(1-x)^2+(V*G)/(R*T)*(surd(1-x,3)-(1-x)/2))*(R*T/V), x);

plot(psi(x), x=0..3, title="The plot of function y=psi(x)");

q:=y->fsolve(psi(x)=y, x):

plot('q'(y), y=-.9e9..-.1.4e9);

If this is important to you, then you can specify  the corresponding function, called P

P:=unapply(n!/(n-r)!, n,r);
 P(10,6);  P(n,r);

                        n!
P := (n, r) -> --------
                     (n - r)!

151200

   n!
--------
(n - r)!

It is easy to see that the product of the members of which are equidistant from the ends is 4, for example

simplify((sqrt(3)+tan(Pi/180))*(sqrt(3)+tan(29*Pi/180)));
                                           4

So, we obtain

mul(simplify((sqrt(3)+tan(Pi*k/180))*(sqrt(3)+tan((30-k)*Pi/180))), k=1..14)*simplify(sqrt(3)+tan(15*Pi/180));
                                     536870912

z:=(x,y)->5/(1+x^2+y^2):

A:=plot3d(z(x,y), x=-4..7, y=-4..4, numpoints=3000):

x:=t->2+3*cos(t): y:=t->-1+2*sin(t):

N:=100:

B:=seq(plots[spacecurve]([x(t),y(t),z(x(t), y(t))], t=0..2*Pi*k/N, color=black, thickness=3), k=0..N):

bug:=seq(plottools[sphere]([x(2*Pi*k/N),y(2*Pi*k/N),z(x(2*Pi*k/N), y(2*Pi*k/N))], 0.15, style=surface, color=red),  k=0..N):

C:=seq(plots[display](A, B[k], bug[k]), k=1..N+1):

plots[display](C, insequence=true, axes=normal, view=[-4.7..7.7, -4.7..4.7, -1..5.7]);

Your task is similar to the problem of wandering drunk sailor. Do a keyword search drunken sailor.

Vector(10,[seq(sin(i), i=-0.5..-0.1, 0.1), seq(sin(i), i=0.1..0.5, 0.1)]);

or

Vector(10, [seq(sin(i), i in {seq(-0.5..0.5, 0.1)} minus {0.})]);

subs(a=`3`, sqrt(a)/a);

To find the center of rotation is sufficient to use two points:

restart; with(geometry):

point(X, x0, y0), point(A, 4, -2), point(B, 5, -4):   #  X - the center of rotation

L:=coordinates(rotation(A1, A, alpha, 'counterclockwise', X)):

M:=coordinates(rotation(B1, B, alpha, 'counterclockwise', X)):

solve({L[1]=4, L[2]=2, M[1]=6, M[2]=3});

You can use  seq  command.

Example:

X:=[seq(i, i=0..20)]:
Y:=[seq((-1)^i*x, i=0..19)]:
y:=piecewise(seq(op([X[i]<x and x<=X[i+1], Y[i]]), i=1..20)):
plot(y, x=0..20);

The last line should be

dsolve({eq, bcs}, v(x));

Matrix([seq([seq(a[k], k = x+m .. y+m)], m = 0 .. 2)]);

a:=1:  c:=2:

fsolve(5*b^5+(60-5*a)*b^4+(125+50*c-80*a)*b^3+(594*c-445*a-775)*b^2+(2324*c-1005*a-3270)*b+3000*c-750*a-3000=0, b);

                                          -5.950891678, -4.378215842, -3.371729567

PS. We find the real roots. If also complex roots are needed, then

fsolve(5*b^5+(60-5*a)*b^4+(125+50*c-80*a)*b^3+(594*c-445*a-775)*b^2+(2324*c-1005*a-3270)*b+3000*c-750*a-3000=0, b, complex); 

                                          -5.950891678, -4.378215842, -3.371729567,
                              1.350418543 - 1.816274549 I, 1.350418543 + 1.816274549 I

Consider the function  f:=x->sqrt(a* x + b) +  sqrt(c*x + d) . If  a*c>=0  then  f   is a monotonic function. Therefore, the equation   sqrt(a* x + b) +  sqrt(c*x + d) = m  can not have two solutions. So the necessary condition is   a*c<0 . For definiteness, let  a>0 . From these conditions and from the conditions a*x+b>=0,  c*x+d>=0  we get  -b/a<=x<=-d/c . Therefore, in the limited ranges of parameters  a, b, c, d  all the solutions can be found by the usual brute force.

The following code finds all equations with integers  a=1 .. 10,  b=-10 .. 10,  c=-10 .. -1,  d=-10 .. 10, each of which has exactly two integer solutions:

L:=[]:

for a to 10 do

for b from -10 to 10 do

for c from -10 to -1 do

for d from -10 to 10 do

if -b/a<=-d/c then s:=floor(-b/a): t:=ceil(-d/c):

M:=[]:

for x from s to t do

u:=sqrt(a*x+b): v:=sqrt(c*x+d):

if type(u, integer) and type(v, integer) then M:=[op(M), [x, u+v]]: fi:

od:

if nops(M)=2 and M[1,2]=M[2,2] then L:=[op(L), [a,b,c,d,M[2,2],[M[1,1],M[2,1]]]]: fi: fi:

od: od: od: od:

nops(L);

L[1..50];

ListTools[Search]([1, 5, -1, 8, 5, [-1,4]], L); 

For example, the list  [1, -9, -1, 10, 1, [9, 10]]  corresponds to the equation  sqrt(x-9)+sqrt(-x+10)=1  with the roots 9 and 10. The list  L  contains all the solutions in the specified ranges. Received 319 solutions. Displayed first 50 solutions. The original equation is also in the list  L  at position 82.

Check your syntax! I have no problem: