Numpy Dot Product of two 2-d arrays in numpy to get 3-d array

One way would be to use np.einsum, which allows you to specify what you want to happen to the indices:

>>> np.einsum('ik,jk->kij', puy, puy2)
array([[[ 0,  0,  0],
        [ 0, 16, 32]],

       [[ 1,  5,  9],
        [ 5, 25, 45]],

       [[ 4, 12, 20],
        [12, 36, 60]],

       [[ 9, 21, 33],
        [21, 49, 77]]])
>>> np.allclose(np.einsum('ik,jk->kij', puy, puy2), zz)
True

Here's another way with broadcasting -

(puy[None,...]*puy2[:,None,:]).T