I am aware about that inverting a matrix to solve a linear system is not a good idea, since it is not as accurate and as efficient as directly solving the system or using LU, Cholesky or QR decomposition. However, I have not been able to check this with a practical example. I have tried this code (in MATLAB) M = 500; A = rand(M,M); A = real(expm(1i*(A+A.'))); b = rand(M,1); x1 = A\b; x2 = inv(A)*b; disp(norm(b-A*x1)) disp(norm(b-A*x2)) and the residuals are always of the same order (10^-13). Could someone provide a practical example in which inv(A)*b is much less inaccurate than A\b?