interpolation with multi-dimensional grid data

Suppose we want to interpolate a multi-dimensional function

y = f(x1, x2, x3...) 

and we have numerical table of values at grid points of each axis. 

This interpolation can be done by using scipy.interpolate.interpn 

However, the data have to be prepared in specific forms. 

scipy.interpolate.interpn(points, values, xi, method='linear', bounds_error=True, fill_value=nan)

(1) "points" are tuple of grid points in each axis. 

     points = ( [array of possible grid points in 1st axis], [array of possible grid points in 2nd axis],...) 

     in other words, points=(x1,x2,...xn) such that  xn[i] gives i-th grid value in n-th axis. 

(2) "values" are n-dimensional array such that 

     value[ i1][i2][i3]...[in] = f( x1[i1],  x2[i2],..., xn[in] ) 

(3) "xi" is a coordinate to sample the value