Thanks a lot, Chris!! This fix indeed increases significantly the performance! However it still much less performant than the Python version and scales really badly with the size of the system. For example, for n = 20 :
julia> @time solve(prob, MethodOfSteps(BS3()), saveat=0.01, reltol=0, abstol=1e-5)
119.585544 seconds (11.80 M allocations: 1.706 GiB, 3.26% gc time)
while in the Python version it takes ~ 0.007 seconds.
My real need is to run it for n = 180 so I’ll need to improve the performance much more.