If you solve the problem for one return value, you can find two functions f_1(x,y) and f_2(x,y) by interpolation and compose your function as f(x, y) = [f_1(x,y), f_2(x,y)] . Just choose any method to solve the interpolation function that suits your problem.
For the actual problem of interpolation in two dimensions, there are many ways to handle this. If this is just what you need, you can go with linear interpolation. If you're fine with piecewise functions, you can use Bezier curves or splines. Or, if the data is homogeneous, you can get away with simple polynomial interpolation (well, not quite trivial when in 2D, but easy enough).
EDIT: additional information and some links.
A piecewise solution is possible using bilinear interpolation (Wikipedia) .
For polynomial interpolation, if your data is in a grid, you can use the following algorithm (I can not find a link for it, this is from memory).
If the data points are on the grid k by l , rewrite your polynomial as follows:
f(x,y) = cx_1(x)*y^(k-1) + cx_2(x)*y^(k-2) + ... + cx_k(x)
Here, each coefficient cx_i(x) also a polynomial of degree l . The first step is to search for polynomials k degree l by interpolating each row or column of the grid. When this is done, you have coefficient coefficients l (or, in other words, polynomials l ) as interpolation points for each polynomial cx_i(x) like cx_i(x0) , cx_i(x1) , ..., cx_i(xl) ( giving you just 1 * k points). Now you can define these polynomials using the above constants as interpolation points that give you the result f(x,y) .
The same method is used for Bezier curves or splines. The only difference is that you use control points instead of polynomial coefficients. First you get a set of splines that will generate your data points, and then you interpolate the control points of these intermediate curves to get the control points of the surface curve.
Let me add an example to clarify the above algorithm. Let have the following data:
0,0 => 1 0,1 => 2 1,0 => 3 1,1 => 4
We start by fitting two polynomials: one for data points (0,0) and (0,1), and the other for (1, 0) and (1, 1):
f_0(x) = x + 1 f_1(x) = x + 3
Now we interpolate in the other direction to determine the coefficients. When we read these polynomial coefficients vertically, we need two polynomials. One evaluates to 1 at 0 and 1; and another that evaluates to 1 at 0 and 3 at 1:
cy_1(y) = 1 cy_2(y) = 2*y + 1
If we combine them in f(x,y) , we get:
f(x,y) = cy_1(y)*x + cy_2(y) = 1*x + (2*y + 1)*1 = x + 2*y + 1