Does everyone have a good idea of ​​the work of the Draghilev method? For example, in this answer https://www.mapleprimes.com/questions/235407-The-Second-Example-Of-Finding-All-Solutions#answer291268
( https://www.mapleprimes.com/questions/235407-The-Second-Example-Of-Finding-All-Solutions )
there was a very successful attempt by Rouben Rostamian to calculate the line of intersection of surfaces without applying the Draghilev method.
Let now not 3d, but 8d. And how will the solve command work in this case? Imagine that aij ((i=1..7,j=1..8)) are partial derivatives, and xj (,j=1..8) are derivatives, as in the above example. f8 is responsible for the parametrization condition.
 

restart;
 f1 := a11*x1+a12*x2+a13*x3+a14*x4+a15*x5+a16*x6+a17*x7+a18*x8; 
 f2 := a21*x1+a22*x2+a23*x3+a24*x4+a25*x5+a26*x6+a27*x7+a28*x8; 
 f3 := a31*x1+a32*x2+a33*x3+a34*x4+a35*x5+a36*x6+a37*x7+a38*x8; 
 f4 := a41*x1+a42*x2+a43*x3+a44*x4+a45*x5+a46*x6+a47*x7+a48*x8;
 f5 := a51*x1+a52*x2+a53*x3+a54*x4+a55*x5+a56*x6+a57*x7+a58*x8; 
 f6 := a61*x1+a62*x2+a63*x3+a64*x4+a65*x5+a66*x6+a67*x7+a68*x8; 
 f7 := a71*x1+a72*x2+a73*x3+a74*x4+a75*x5+a76*x6+a77*x7+a78*x8;
 f8 := x1^2+x2^2+x3^2+x4^2+x5^2+x6^2+x7^2+x8^2-1; 
allvalues(solve({f1, f2, f3, f4, f5, f6, f7, f8}, {x1, x2, x3, x4, x5, x6, x7, x8}))

And this is how the Draghilev method works in this case.
 

restart; with(LinearAlgebra):
f1 := a11*x1+a12*x2+a13*x3+a14*x4+a15*x5+a16*x6+a17*x7+a18*x8; 
f2 := a21*x1+a22*x2+a23*x3+a24*x4+a25*x5+a26*x6+a27*x7+a28*x8; 
f3 := a31*x1+a32*x2+a33*x3+a34*x4+a35*x5+a36*x6+a37*x7+a38*x8; 
f4 := a41*x1+a42*x2+a43*x3+a44*x4+a45*x5+a46*x6+a47*x7+a48*x8; 
f5 := a51*x1+a52*x2+a53*x3+a54*x4+a55*x5+a56*x6+a57*x7+a58*x8;
f6 := a61*x1+a62*x2+a63*x3+a64*x4+a65*x5+a66*x6+a67*x7+a68*x8; 
f7 := a71*x1+a72*x2+a73*x3+a74*x4+a75*x5+a76*x6+a77*x7+a78*x8;

n := 7; 
x := seq(eval(cat('x', i)), i = 1 .. n+1):
 F := [seq(eval(cat('f', i)), i = 1 .. n)]: 
A := Matrix(nops(F), nops(F)+1):
 for j to nops(F) do for i to nops(F)+1 do A[j, i] := op(1, op(i, op(j, F))) 
end do:
           end do: 

# b[i] and b[n+1] are solutions of a linear homogeneous system and at the 
# same time they serve as the right-hand sides of an autonomous ODE.
for i to n do

 b[i] := Determinant(DeleteColumn(ColumnOperation(A, [i, n+1]), n+1)) 
                                                                     end do:
 b[n+1] := -Determinant(DeleteColumn(A, n+1)):

Only the original seven linear homogeneous equations with eight variables are needed. We solve them according to Cramer's rule, and in order to have uniqueness when solving the ODE, we use a point on the curve (according to the theory). (This point is sought in any convenient way.)
If we want to get a parameterization, then additionally, directly in dsolve, we can add the following:
 

for i to n do 
b[i] := simplify(Determinant(DeleteColumn(ColumnOperation(A, [i, n+1]), n+1))) 
end do:
b[n+1] := simplify(-Determinant(DeleteColumn(A, n+1))); 
deqs := seq(diff(x[i](s), s) = b[i]/(b[1]^2+b[2]^2+b[3]^2+b[4]^2+b[5]^2+b[6]^2+b[7]^2+b[8]^2)^.5, i = 1 .. n+1):