Inversion and multiplication of sparse matrices

Coming from Matlab, I was expecting the same behavior too.

As far as I can tell, Julia’s \ accepts as right-hand side b only vectors like rand(4) or single-column matrices like rand(4,1). The result is of the same type as b.
The vector form seems to be preferred, not all linear solvers accept a single-column matrix.

See also Avoiding memory allocation for solves with multiple right-hand sides - General Usage / Performance - JuliaLang

using LinearAlgebra, SparseArrays

function solve_multiRHS!(A, b)
    lu!(A)
    for i = 1:length(b)
        b[i] =  A\b[i]
    end
    return b
end

A = rand(4,4)
b = Vector{Float64}[]
for i = 1:3
    push!(b, rand(4))
end 

solve_multiRHS!(A, b)

The sparse operator, when needed, is then applied to A and the individual vectors of b.