This code does exactly what I want it to, however I want to try and get rid of the nested loop to make it more pythonic. I have been trying to somehow use broadcasting, including playing with np.newaxis, but cannot produce the same result.

    M1 = np.array([[11,11,11],[12,12,12],[13,13,13]])
    M2 = np.array([[21,21,21],[22,22,22],[23,23,23],[24,24,24]])
    
    m1_rows = M1.shape[0]
    m2_rows = M2.shape[0]
    
    d = np.empty((m1_rows,m2_rows))
    
    for i in range(m1_rows):
        for j in range(m2_rows):
            d[i,j] = fun(M1[i],M2[j])

Some additional details:

M1 and M2 will always be 2 dimensional numpy arrays. They will always have the same number of columns but rows can vary.

fun() takes a row from each matrix and returns a number. Similar to np.dot()