(%i2) eq : 3 * x - y - z + t = 2 $
sol : linsolve([ eq],[ x, y, z, t]) ;
(sol) [ x = %r3 %r2 + %r1 2 3 , y = %r3 , z = %r2 , t = %r1 ]

Expresemos las soluciones en forma matricial:
Un punto de la variedad lineal solución del sistema sería el dado con los parámetros igual a cero, y los vectores:

(%i6) P : ev([ x, y, z, t], ev( sol, %rnum_list[ 1] = 0, %rnum_list[ 2] = 0, %rnum_list[ 3] = 0)) $
v1 : ev([ x, y, z, t], ev( sol, %rnum_list[ 1] = 1, %rnum_list[ 2] = 0, %rnum_list[ 3] = 0)) - P $
v2 : ev([ x, y, z, t], ev( sol, %rnum_list[ 1] = 0, %rnum_list[ 2] = 1, %rnum_list[ 3] = 0)) - P $
v3 : ev([ x, y, z, t], ev( sol, %rnum_list[ 1] = 0, %rnum_list[ 2] = 0, %rnum_list[ 3] = 1)) - P $

La forma matricial es:

(%i7) print( matrix([ x],[ y],[ z],[ t]), "=", transpose( matrix( P)), "+ λ",
   transpose( matrix( v1)), "+ μ",
   transpose( matrix( v2)), "+ ν",
   transpose( matrix( v3))) $ ;
[ x y z t ] = [ 2 3 0 0 0 ] + λ [ 1 3 0 0 1 ] + μ [ 1 3 0 1 0 ] + ν [ 1 3 1 0 0 ]

El vector (a,3,2,1) pertenece al subespacio director de la variedad si es combinación lineal de los vectores de una base. Lo que es equivalente a que el rango de la matriz rank(matrix([a,3,2,-1],v1,v2,v3))=3.

Para que el rango sea 3, necesitamos que el determinante de orden 4 sea cero:

(%i8) solve( determinant( matrix([ a, 3, 2, - 1], v1, v2, v3)), a) ;
(%o8) [ a = 2 ]

Created with wxMaxima.