solveで,Point3Dが使えました。
4点 A（1，1，3），B（－1，1，－1），C（1，－4，－2），D（－2，－1，1）

公式サイト

<問題>数学（Z）問1

<解答例>

上と同じです。大学入試数学問題集成＞【1】

Pycharmで

答えは、あっていると思います。

from sympy import *
var('s t u')
def myInnerProduct(myPA,myPB):
    myA=myPoint3DtoMatrix(myPA)
    myB=myPoint3DtoMatrix(myPB)
    mySum=0
    for i in range(0, shape(myA)[1]):
        mySum = mySum+myA[i]*myB[i]
    return mySum
def myPoint3DtoMatrix(myPoint3D):
    myx = Matrix(1, 3, range(3))
    myx[0]=myPoint3D.x
    myx[1]=myPoint3D.y
    myx[2]=myPoint3D.z
    return myx
myZahyo=((0,0,0),(1,1,3),(-1,1,-1),(1,-4,-2),(-2,-1,1))
O=Point3D(myZahyo[0])
A=Point3D(myZahyo[1])
B=Point3D(myZahyo[2])
C=Point3D(myZahyo[3])
D=Point3D(myZahyo[4])
ans=solve(O-D-(s*(A-D)+t*(B-D)+u*(C-D)),[s,t,u])
print("#〔1〕s,t,u:",ans[s],ans[t],ans[u])
#
var('v w')
P=A+v*(B-A)
Q=C+w*(D-C)
ans=solve([myInnerProduct(B-A,Q-P),myInnerProduct(Q-P,D-C)],[v,w])
mySubs={v:ans[v],w:ans[w]}
myFutosiki= 0 <= v.subs(mySubs) <= 1
print("#〔2〕線分上:",type(myFutosiki),myFutosiki)
print("#    PQ   :",P.distance(Q).subs(mySubs),
               "P:",P.x.subs(mySubs),P.y.subs(mySubs),P.z.subs(mySubs),
               "Q:",Q.x.subs(mySubs),Q.y.subs(mySubs),Q.z.subs(mySubs))
#
var('x')
R=C+x*(D-C)
cosARB=myInnerProduct(A-B,R-B)/(Point(B).distance(A)*Point(B).distance(R))
ans=solve(diff(cosARB))
mySubs={x:ans[0]}
print("#〔3〕ABR  :",ans[0],cosARB.subs({x:ans[0]}),
               "R:",R.x.subs(mySubs),R.y.subs(mySubs),R.z.subs(mySubs))
#〔1〕s,t,u: 3/10 1/2 1/5
#〔2〕線分上: <class 'sympy.logic.boolalg.BooleanTrue'> True
#    PQ   : 5*sqrt(14)/7 P: -3/7 1 1/7 Q: -13/7 -8/7 6/7
#〔3〕ABR  : 5/4 sqrt(39)/13 R: -11/4 -1/4 7/4

参考

wolframalphaでグラフを作図

x = 5/4 のとき, max{(3 sqrt(5) x)/(5 sqrt((1 - 3 x)^2 + (5 - 3 x)^2 + (3 x - 2)^2))} = sqrt(3/13)

FreeCADのマクロで

上のsympy計算に、FreeCAD作図を追加しました。
操作で、線2本を、赤色と緑色にしました。

import FreeCAD
import Part
import Draft
import Mesh
#########################################################################################################
# 計算
from sympy import *
var('s t u')
def myInnerProduct(myPA,myPB):
    myA=myPoint3DtoMatrix(myPA)
    myB=myPoint3DtoMatrix(myPB)
    mySum=0
    for i in range(0, shape(myA)[1]):
        mySum = mySum+myA[i]*myB[i]
    return mySum
def myPoint3DtoMatrix(myPoint3D):
    myx = Matrix(1, 3, range(3))
    myx[0]=myPoint3D.x
    myx[1]=myPoint3D.y
    myx[2]=myPoint3D.z
    return myx
myZahyo=((0,0,0),(1,1,3),(-1,1,-1),(1,-4,-2),(-2,-1,1))
O=Point3D(myZahyo[0])
A=Point3D(myZahyo[1])
B=Point3D(myZahyo[2])
C=Point3D(myZahyo[3])
D=Point3D(myZahyo[4])
ans=solve(O-D-(s*(A-D)+t*(B-D)+u*(C-D)),[s,t,u])
print("#〔1〕s,t,u:",ans[s],ans[t],ans[u])
#
var('v w')
P=A+v*(B-A)
Q=C+w*(D-C)
ans=solve([myInnerProduct(B-A,Q-P),myInnerProduct(Q-P,D-C)],[v,w])
mySubs={v:ans[v],w:ans[w]}
P=Point3D(P.x.subs(mySubs),P.y.subs(mySubs),P.z.subs(mySubs))
Q=Point3D(Q.x.subs(mySubs),Q.y.subs(mySubs),Q.z.subs(mySubs))
myFutosiki= 0 <= v.subs(mySubs) <= 1
print("#〔2〕線分上:",type(myFutosiki),myFutosiki)
print("#     PQ    :",P.distance(Q),"P:",P,"Q:",Q)
#
var('x')
R=C+x*(D-C)
cosARB=myInnerProduct(A-B,R-B)/(Point(B).distance(A)*Point(B).distance(R))
ans=solve(diff(cosARB))
mySubs={x:ans[0]}
R=Point3D(R.x.subs(mySubs),R.y.subs(mySubs),R.z.subs(mySubs))
print("#〔3〕ABR    :",ans[0],cosARB.subs({x:ans[0]}),
                  "R:",R)
#〔1〕s,t,u: 3/10 1/2 1/5
#〔2〕線分上: <class 'sympy.logic.boolalg.BooleanTrue'> True
#    PQ   : 5*sqrt(14)/7 P: -3/7 1 1/7 Q: -13/7 -8/7 6/7
#〔3〕ABR  : 5/4 sqrt(39)/13 R: -11/4 -1/4 7/4
# #########################################################################################################
# 3D作図
def myXYZ2Txt(A):
    return '(' + str(A.x) + ',' + str(A.y) + ',' + str(A.z) + ')'
def myTxt_XYZ(A,myWedgei):
    P5x=float(A.x)
    P5y=float(A.y)
    P5z=float(A.z)
    p5 = FreeCAD.Vector(P5x, P5y, P5z)
    myText = Draft.makeText(myWedgei, p5)
    myText.Label = myWedgei
    FreeCADGui.ActiveDocument.ActiveObject.FontSize = '0.4mm'
def myLine(A,B):
    Ax,Ay,Az=float(A.x),float(A.y),float(A.z)
    Bx,By,Bz=float(B.x),float(B.y),float(B.z)
    pl = FreeCAD.Placement()
    pl.Rotation.Q = (0.4247081540122249, 0.17592004639554645, 0.33985110062924484, 0.8204732460821097)
    pl.Base = FreeCAD.Vector(-3.9166066876399563, -2.1670824762243774, 1.7495260956243028)
    points = [FreeCAD.Vector(Ax,Ay,Az), FreeCAD.Vector(Bx,By,Bz)]
    line = Draft.make_wire(points, placement=pl, closed=False, face=True, support=None)
    Draft.autogroup(line)
myTxt_XYZ(O,"O"+myXYZ2Txt(O))
myTxt_XYZ(A,"A"+myXYZ2Txt(A))
myTxt_XYZ(B,"B"+myXYZ2Txt(B))
myTxt_XYZ(C,"C"+myXYZ2Txt(C))
myTxt_XYZ(D,"D"+myXYZ2Txt(D))
myTxt_XYZ(P,"P"+myXYZ2Txt(P))
myTxt_XYZ(Q,"Q"+myXYZ2Txt(Q))
myTxt_XYZ(R,"R"+myXYZ2Txt(R))
myLine(A,B)
myLine(C,D)
myLine(P,Q)
myLine(B,R)
myLine(D,R)
doc = App.activeDocument()
App.ActiveDocument.addObject("App::Origin", "Origin")
App.ActiveDocument.getObject('Origin').Visibility = True
App.ActiveDocument.recompute()
Gui.activeDocument().activeView().viewAxonometric()
Gui.SendMsgToActiveView("ViewFit")