A = matrix([[1, 3, 3], [2, 6, 9], [-1, -3, 3]])
print A
b = vector([1, 5, 5])
x = A.solve_right(b)
print x
A * x
x, y, z = var('x, y, z')
solve([x+3*y+3*z==1, 
       2*x+6*y+9*z==5,
       -x-3*y+3*z==5], x, y, z)