I have a function fun(vec(1:n),Nmax) with two outputs: (i) a matrix g(1:n,1:Nmax) and (ii) a vector tau(1:Nmax), as Nmax eigenfunctions and eigenvalues of a matrix constructed within fun. Now I want to get this output for all values withing a vector vec2(1:n). The simplest way is a for-loop
[g(1:n,1:Nmax,1:n),tau(1:n,1:Nmax)] = fun(vec(1:n),Nmax,vec2(i))
This is however slow. Replacing it with a parfor-loop
[g(:,:,i),tau(i,:)] = fun(vec(:),Nmax,vec2(i))
is quite a bit faster. I was wondering if there is a way to accelerate this even further by threading over the vector vec2? Somehow I cannot find the correct way without rewriting the function fun in several ways. Maybe the problem is that g is already a matrix, therefore fun(vec,Nmax,vec2') does not evaluate things in the correct dimension?