Convert solution set to vector notation?

edit

I have the following matrix and RREF:

[
[ 2  1  0 -1  0| 4]    [   1    0    0   -1 -1/2|   0]
[ 0  1  0  1  1| 4]    [   0    1    0    1    1|   4]
[ 1  0 -1  2  0| 0],   [   0    0    1   -3 -1/2|   0]
]

Is it possible to convert the set of solutions:

  x == 1/2*u + w
  y == -u - w + 4
  z == 1/2*u + 3*w

into vector form:

[ x ]   [ 0 ]   [  1 ]      [ 1/2 ]   |   
[ y ]   [ 4 ]   [ -1 ]      [  -1 ]   |
[ z ] = [ 0 ] + [  3 ] w +  [ 1/2 ]u  |   w, u ∈ ℝ
[ w ]   [ 0 ]   [  1 ]      [   0 ]   |
[ u ]   [ 0 ]   [  0 ]      [   1 ]   |

Is this possible with Sagemath?

The solutions are created using:

def my_solve(augmented_matrix, vars=None):

    A = augmented_matrix[:, :-1]
    Y = augmented_matrix[:, -1]

    m, n = A.dimensions()
    p, q = Y.dimensions()

    if m!=p:
        raise RuntimeError("The matrices have different numbers of rows")

    if vars and len(vars) != n:
        raise RuntimeError(f"Provided variables '{vars}' != number of columns '{n}'")

    if vars:
        # don't include the free variables in solve
        X = vector([var(vars[i]) for i in [0..n-1]])
        X_pivots = vector([var(X[i]) for i in [0..n-1] if i in A.pivots()])
        X_free = vector([var(X[i]) for i in [0..n-1] if i not in A.pivots()])
    else:
        X = vector([var(f"x_{i}") for i in [0..n-1]])
        X_pivots = vector([var(f"x_{i}") for i in [0..n-1] if i in A.pivots()])
        X_free = vector([var(f"x_{i}") for i in [0..n-1] if i not in A.pivots()])

    sols = []
    for j in range(q):
        system = [A[i]*X==Y[i,j] for i in range(m)]
        sols += solve(system, *X_pivots)

    return sols, X, X_pivots, X_free

M = matrix(QQ, 3, [2,1,0,-1,0,   0,1,0,1,1,   1,0,-1,2,0])
v = vector(QQ, [4,4,0])
Maug = M.augment(v, subdivide=True)
my_solve(Maug, var('x y z w u'))[0]