If it’s stable with ode45, then it should be stable with Tsit5, DP5, etc. While your style of coding doesn’t have to change, I can’t say it’s the easiest style to debug since the indexing is a little… intense.
Your initial condition is very constant, so that’s probably not a good place to be checking the correctness of the derivative evaluation, since applying some of those operators might give be giving zeros due to the constant-ness. I would generate a random vector and use the same random initial condition in Julia and MATLAB and then invoke the f call and see if the derivatives are the same. My guess is that they aren’t.
(Note: this coding style does heavily allocate though, so the speed will not be top notch because of all of those extra array creations. You could definitely implement this derivative function in a different way that would speed it up quite a bit)