Speed up Numpy Meshgrid Command

def meshgrid(x,y):
    """
    Return coordinate matrices from two coordinate vectors.

    Parameters
    ----------
    x, y : ndarray
        Two 1-D arrays representing the x and y coordinates of a grid.

    Returns
    -------
    X, Y : ndarray
        For vectors `x`, `y` with lengths ``Nx=len(x)`` and ``Ny=len(y)``,
        return `X`, `Y` where `X` and `Y` are ``(Ny, Nx)`` shaped arrays
        with the elements of `x` and y repeated to fill the matrix along
        the first dimension for `x`, the second for `y`.

    See Also
    --------
    index_tricks.mgrid : Construct a multi-dimensional "meshgrid"
                         using indexing notation.
    index_tricks.ogrid : Construct an open multi-dimensional "meshgrid"
                         using indexing notation.

    Examples
    --------
    >>> X, Y = np.meshgrid([1,2,3], [4,5,6,7])
    >>> X
    array([[1, 2, 3],
           [1, 2, 3],
           [1, 2, 3],
           [1, 2, 3]])
    >>> Y
    array([[4, 4, 4],
           [5, 5, 5],
           [6, 6, 6],
           [7, 7, 7]])

    `meshgrid` is very useful to evaluate functions on a grid.

    >>> x = np.arange(-5, 5, 0.1)
    >>> y = np.arange(-5, 5, 0.1)
    >>> xx, yy = np.meshgrid(x, y)
    >>> z = np.sin(xx**2+yy**2)/(xx**2+yy**2)

    """
    x = asarray(x)
    y = asarray(y)
    numRows, numCols = len(y), len(x)  # yes, reversed
    x = x.reshape(1,numCols)
    X = x.repeat(numRows, axis=0)

    y = y.reshape(numRows,1)
    Y = y.repeat(numCols, axis=1)
    return X, Y

def meshgrid2(*arrs):
    arrs = tuple(reversed(arrs))  #edit
    lens = map(len, arrs)
    dim = len(arrs)

    sz = 1
    for s in lens:
        sz*=s

    ans = []    
    for i, arr in enumerate(arrs):
        slc = [1]*dim
        slc[i] = lens[i]
        arr2 = asarray(arr).reshape(slc)
        for j, sz in enumerate(lens):
            if j!=i:
                arr2 = arr2.repeat(sz, axis=j) 
        ans.append(arr2)

    return tuple(ans)

>>> xi, yi, zi = [np.arange(3) for i in range(3)]
>>> xx, yy, zz = np.broadcast_arrays(xi,yi[:,np.newaxis],zi[:,np.newaxis,np.newaxis])
>>> xx.shape
(3, 3, 3)
>>> xx.strides
(0, 0, 8)

def ndmesh(*args):
   args = map(np.asarray,args)
   return np.broadcast_arrays(*[x[(slice(None),)+(None,)*i] for i, x in enumerate(args)])

X, Y = np.mgrid[xmin:xmax:100j, ymin:ymax:100j]
positions = np.vstack([X.ravel(), Y.ravel()])

g = meshgrid2(x, y, z)
positions = np.vstack(map(np.ravel, g))

In [23]: import numpy as np

In [34]: from numpy import asarray

In [35]: x = np.random.rand(100,1)

In [36]: y = np.random.rand(100,1)

In [37]: z = np.random.rand(100,1)

In [38]: def meshgrid2(*arrs):
   ....:     arrs = tuple(reversed(arrs))
   ....:     lens = map(len, arrs)
   ....:     dim = len(arrs)
   ....:     sz = 1
   ....:     for s in lens:
   ....:        sz *= s
   ....:     ans = []
   ....:     for i, arr in enumerate(arrs):
   ....:         slc = [1]*dim
   ....:         slc[i] = lens[i]
   ....:         arr2 = asarray(arr).reshape(slc)
   ....:         for j, sz in enumerate(lens):
   ....:             if j != i:
   ....:                 arr2 = arr2.repeat(sz, axis=j)
   ....:         ans.append(arr2)
   ....:     return tuple(ans)

In [39]: g = meshgrid2(x, y, z)

In [40]: %timeit pos = np.vstack(map(np.ravel, g)).T
100 loops, best of 3: 7.26 ms per loop

In [41]: %timeit zip(*(x.flat for x in g))
1 loops, best of 3: 264 ms per loop

from numpy import array, meshgrid

def mesh(foo):
  return array(meshgrid(*foo)).T.reshape(-1, foo.shape[0])

print(mesh(array(([1, 2, 3], [4, 5]))))
print(mesh(array(([1, 2, 3], [4, 5], [6, 7]))))