Let's say you know x1, x2 and y1=y(x1), y2=y(x2) and gradients dy/dx at x1 and dy/dx at x2. That's 4 knowns. Thus, you can uniquely define the parameters in a cubic curve y=a+bx+cx^2+dx^3, which has 4 unknowns, that interpolate the function between x1 and x2. The interpolated cubic spline can be made to give the correct y values and derivatives at x1, x2. OK, for 2 dimensions- each point has 1 y value and it also has a gradient, where the gradient now has two components. Thus, each point has 3 knowns attached to it. If you have 4 points- say on the corners of a square, where each point has known y values and known gradients, then there are 12 knowns altogether. Unfortunately, a cubic equation in two dimensions has (the triangular number) 10 coefficients. So there's a mismatch between the number of knowns and unknowns. One solution seems to be to use a quartic interpolation, which has (the triangular number) 15 coefficients. Then you could use a 5 points, (for example, corners of a square, with a point in the middle) with 5*3 = 15 knowns, and there would be an even match between knowns and unknowns. 5 points in 3 dimensions also seems to work. Is there a general way of solving this problem for N dimensions? I'm sure this stuff must have been worked out.