Comment by geysersam

You're right but if you want `y` to have the same shape as `x` I think you also need to slice it after

    y = linalg.solve(A, x[..., None])[..., 0]

I don't really mind numpy syntax much, it gets the job done in most scenarios. Numba complements it really well when the code is easier to express like a couple of nested loops.

I used Julia for a while, but found it easier to predict the performance of python+numpy+numba, Julia has a few footguns and the python ecosystem is just insanely polished.