Matrix solution with Numpy: no solution?
If you take pen and paper (or use numpy.linalg.matrix_rank) and calculate the ranks of the coefficient and the augmented matrix you'll see that, according to Rouché–Capelli theorem, there is no solution. So Wolfram is right.
NumPy searches for numerical solution of equation using LU decomposition. Without going into details LU decomposition involves divisions, divisions in FP arithmetic can lead to significant errors.
If you check:
a = np.matrix('1 4 1; 4 13 7; 7 22 13')b = np.matrix('0;0;1')c = np.linalg.solve(a,b)np.linalg.norm(a * c - b)## 2.0039024427351748you'll see that NumPy solution is far from being correct.