Before using PETSc, I tried the simple way: sol = A^-1 * rhs, where sol and rhs are 2D arrays. Then I extract each solution as sol(:, i). However, the solution obtained this way are completely wrong unless rhs has only 1 column. For example, one solution looks like this:
