(参考)
sympyで
fullscript.py
from sympy import *
var("a11 a12 a13 a21 a22 a23 a31 a32 a33 b1 b2 b3")
var("x y z ")
var("A B C D")
print(solve([2*x+5*y+7*z-4,4*x+13*y+20*z-11,8*x+29*y+50*z-29 ]))
A =Matrix(([a11,a12,a13], [a21,a22,a23],[a31,a32,a33]))
B =Matrix(([b1],[b2],[b3]))
C=A.inv() * B
D=C.subs({a11:2,a12:5,a13:7, a21:4, a22:13, a23:20, a31:8,a32:29,a33:50,b1:4, b2:11, b3:29})
print(C)
print(D)
print("x=",D[0],"y=",D[1],"z=",D[2])
# {x: 1, y: -1, z: 1}
# Matrix([[-b3*(-a12*(a11*a23 - a13*a21) + a13*(a11*a22 - a12*a21))/((a11*a22 - a12*a21)*(a11*a33 - a13*a31) - (a11*a23 - a13*a21)*(a11*a32 - a12*a31)) + b1*(a11*a22*((a11*a22 - a12*a21)*(a11*a33 - a13*a31) - (a11*a23 - a13*a21)*(a11*a32 - a12*a31)) - (-a12*(a11*a23 - a13*a21) + a13*(a11*a22 - a12*a21))*(a21*(a11*a32 - a12*a31) - a31*(a11*a22 - a12*a21)))/(a11*(a11*a22 - a12*a21)*((a11*a22 - a12*a21)*(a11*a33 - a13*a31) - (a11*a23 - a13*a21)*(a11*a32 - a12*a31))) + b2*(-a11*a12*((a11*a22 - a12*a21)*(a11*a33 - a13*a31) - (a11*a23 - a13*a21)*(a11*a32 - a12*a31)) + a11*(a11*a32 - a12*a31)*(-a12*(a11*a23 - a13*a21) + a13*(a11*a22 - a12*a21)))/(a11*(a11*a22 - a12*a21)*((a11*a22 - a12*a21)*(a11*a33 - a13*a31) - (a11*a23 - a13*a21)*(a11*a32 - a12*a31)))], [-a11*b3*(a11*a23 - a13*a21)/((a11*a22 - a12*a21)*(a11*a33 - a13*a31) - (a11*a23 - a13*a21)*(a11*a32 - a12*a31)) + b1*(-a21*((a11*a22 - a12*a21)*(a11*a33 - a13*a31) - (a11*a23 - a13*a21)*(a11*a32 - a12*a31)) - (a11*a23 - a13*a21)*(a21*(a11*a32 - a12*a31) - a31*(a11*a22 - a12*a21)))/((a11*a22 - a12*a21)*((a11*a22 - a12*a21)*(a11*a33 - a13*a31) - (a11*a23 - a13*a21)*(a11*a32 - a12*a31))) + b2*(a11*(a11*a23 - a13*a21)*(a11*a32 - a12*a31) + a11*((a11*a22 - a12*a21)*(a11*a33 - a13*a31) - (a11*a23 - a13*a21)*(a11*a32 - a12*a31)))/((a11*a22 - a12*a21)*((a11*a22 - a12*a21)*(a11*a33 - a13*a31) - (a11*a23 - a13*a21)*(a11*a32 - a12*a31)))], [-a11*b2*(a11*a32 - a12*a31)/((a11*a22 - a12*a21)*(a11*a33 - a13*a31) - (a11*a23 - a13*a21)*(a11*a32 - a12*a31)) + a11*b3*(a11*a22 - a12*a21)/((a11*a22 - a12*a21)*(a11*a33 - a13*a31) - (a11*a23 - a13*a21)*(a11*a32 - a12*a31)) + b1*(a21*(a11*a32 - a12*a31) - a31*(a11*a22 - a12*a21))/((a11*a22 - a12*a21)*(a11*a33 - a13*a31) - (a11*a23 - a13*a21)*(a11*a32 - a12*a31))]])
# Matrix([[1], [-1], [1]])
# x= 1 y= -1 z= 1
wolframalphaで
2*x+5*y+7*z-4,4*x+13*y+20*z-11,8*x+29*y+50*z-29
x = 1, y = -1, z = 1