Unfortunately, I get the feeling that we get the same problems in Julia as in C++ ;-(
It is obvious that things need to be done inplace: [1] There is no solution so far and I don't think that kind of problems can be solved easily, except passing "memory location" to functions and mutating them. [2] Template metaprogramming is already a mess in C++. Julia starts with better meta-programming capabilities but I prefer to keep simple things simple (KISS). We mainly need to define : add_linear_combination!(y, lambda1, y1, lambda2, y2, ... , lamdan, yn) and for various values of n which is where meta-programming can be used. I'll implement that in my code. The problem is that `y = y + delta_t * dy_dt` currently makes two copies, > as you have observed. There is an open discussion about in-place assignment > operators [1], which would help. There is also Devectorize.jl [2]. However, > I agree that having something like add!(y, x, beta) would be nice and I > think this also came up in some discussion. Funnily enough we just recently > discussed this strategy for ODE.jl [3] ... >
