WIP: Use rref() in solve_linear_system()

This partially fixes issue 1850. What remains to do in the long term is to add early exit capability to `rref()` (most likely via a flag, because we want to have the possibility to get row-echelon matrix). This is out of scope of this PR. WIP is because there one test failure that is, well, a little confusing.