When writing differential equation solvers, a lot of cases are far too small for vectorization to help and lots of loops are required. Especially adaptive stochastic differential equation solvers which require using data structures (priority queues, stacks) directly in the middle of tight loops, and you have to control heap-allocations to keep it fast. I couldn't find a way how to do this in MATLAB or Python well (I mean, heap-allocated objects tend to ruin this) but Julia was a great fit.