How to apply a function to a Cartesian product efficiently in Python/Numpy?

I have a function that I would like to apply on each element of a cartesian product of a linear space. I do know how to do it with functions of one variable that is defined using lambda, namely using map. Here is an example:

import numpy as np
xpts = np.linspace(0, 1, 5)
fun = lambda p: p**2
arr = np.array(list(map(fun , xpts)))

But with a multivariate function I did not manage to use the map function. Here is an example of what I am doing instead, which is slow:

def fun(x,y):
    return 2*x+y

xpts = np.linspace(0, 1, 5)

# make dict of indexes
count=0
dic=dict()
for j in xpts:
    dic[j]=count
    count+=1

# preallocate array
arr = np.empty([len(xpts)]*2)

for tup in itertools.product(xpts, xpts):

    ind1 = dic[tup[0]]
    ind2 = dic[tup[1]]

    val1 = tup[0]
    val2 = tup[1]

    arr[ind1, ind2]   = fun(val1, val2)

As my function is complicated and the space is large, I am looking for an efficient/scalable way.