Hello, I have 2*N arrays (let's say 2 arrays of length 4) that I combine using np.meshgrid xxA, yyA = np.meshgrid(xA, yA) xxB, yyB = np.meshgrid(xB, yB) which gives me two meshes xx.shape = yy.shape = (4,4) wish represent N-dimensional mesh with 16 elements. Now I want to evaluate a function f on every possible pair of N-dimensional points in the grid, resulting in a 16 x 16 matrix: in a flattened notation, pA = (xxA, yyA) f(pA[1]-pB[1]) f(pA[1]-pB[2]) f(pA[1]-pB[3]) ... f(pA[2]-pB[1]) f(pA[2]-pB[2]) f(pA[2]-pB[3]) ... f(pA[3]-pB[1]) f(pA[3]-pB[2]) f(pA[3]-pB[3]) ... . . . Let's say xA = yA = [1,2,3] and xB = yB = [10,20,30] that gives me a mesh A: (1,3) (2,3) (3,3) (1,2) (2,2) (3,2) (1,1) (2,1) (3,1) and a mesh B alike. My result matrix now should be of size 9 x 9: f( (1,3), (10,30) ) f( (2,3), (20,30) ) f( (3,3), (30, 30) ) f( (1,2), (10,20) ) f( (2,2), (20,20) ) f( (3,2), (30, 20) ) ... f always takes two N-dimensional vectors and returns a scalar. I hope I was able to explain what I want to achieve. What is the best way to do that in numpy? Thanks, Florian