A=matrix([[3, 0, 3, -5], [1, -1, 1, -1], [-2, -1, 4, -2], [2, 1, -1, -1]])
b=vector([-8, -2, 0, -3])
x=A.solve_right(b)
k=A.right_kernel().matrix() # basis of solutions in rows
nrows = k.nrows() # counts the kernel dimension
if nrows > 0: # non unique solution
    t=vector( var("t", n=nrows)) # row of parameters $t_0,\dots$
    show(x + t*k)
else: # unique solution
    show(x)
A.rref() # the reduced EF shows the solutions directly, too