My understanding is that - I admit I haven’t tried this - if you use the ‘bsxfun’ with its ‘fun’ defined as ‘besselj’ and the two arguments “reshape(nVec+1,1,1,[])” and “2*pi*XY”, you will get a 3D array for 'besselForCurrentN' as the result in accordance with nVec+1 varying along the third dimension. An alternative would be to use ‘repmat’ and ‘reshape’ on both nVec+1 and 2*pi*XY in such a way that both become 3D arrays of the same size and then call on ‘besselj’ with these as inputs - I believe ‘besselj’ will accept two 3D arrays if they are the same size.
Whether you can use such a 3D result in a vectorized manner in your further computations is something only you can determine.
