N-dimensional high order polynomial interpolation

Question

I am searching some clues about a complex issue that I am facing, regarding interpolation in a 4D space.

I have a dataset composed by 340 points in a 3-dimensional space (I have three variables - A, B, C - each defined by 340 elements). Each point is identified by a certain value of an output variable. So, generally I have

f(A,B,C) = D

I need to interpolate the dataset in order to predict the value of D for each point in the design space. What I did was to write a small script to obtain the coefficients of the polynomial m through the numpy method linalg.lstsq

def polyfit4d(x,y,z, metric, order):
    ncols = (order + 1)**3
    G = np.zeros((x.size, ncols))
    ij = itertools.product(range(order+1), repeat=3)
    for w, (i,j,k) in enumerate(ij):
        G[:,w] = x**i * y**j * z**k

    m, residuals, rank, s = np.linalg.lstsq(G, metric)

    return m, residuals

Then, I used an evaluating function to obtain the values of the function at all the points of the design space.

def polyval4d(x, y, z, m):
     string = ''
     order = int(math.ceil(((len(m))**(1/3.0))-1))
     ij = itertools.product(range(order+1), repeat=3)
     f = np.zeros_like(x)
     for a, (i,j,k) in zip(m, ij):
     f += a * x**i * y**j * z**k

     return f

Since my design space is 3-dimensional, I passed to the polyval function three 3D matrix with all the points X,Y,Z of the design space. The f is a 3D matrix of outputs D. Each point in this matrix is the value of D calculated evaluating the polynomial, found with polyfit, in each point of the design space (sorry for the tricky sentence).

What I do then is to plot the contour plot of a slice of this 3D design space. I choose one value of Z, and I plot the 2D plane formed by X,Y with the contours levels based on the values of D. The problem is that the result is not what I expect. The contour plot is almost of the same colour, with some variations in one corner.

I searched everywhere on the internet, and also the Python wiki suggests functions that works only for the 2D case. Has anyone ever faced this kind of problem? Do I miss something in the evaluation/definition of this N-dimensional polynomial?

Thanks a lot for your kind attention.

Federico