Help me outperform Matlab in numerical solution of semidiscretized PDE system as much as possible

Thank you for the suggestion. I don’t think that I use matrix multiplication in my DAE code. Even the line of code evaluating the r_R equation is applying element-wise operations on vectors, no?

If this is the case, can I expect a speedup of using Octavian vs. MKL? This might be the case if the solver would be able to make use of this under the hood, I guess.