Interpolation in 2-D co-ordinate system Suppose there is a function f defined on (x, y) such that x, y $\in$ (-$\infty$, $\infty$), but function is not known. Let n data points are given
f($X_0$, $Y_0$) = $Z_0$
f($X_1$, $Y_1$) = $Z_1$
f($X_2$, $Y_2$) = $Z_2$
.......
f($X_n$, $Y_n$) = $Z_n$
These points are not arranged in any particular format, i.e., not equally spaced or any grid formation occurs. Now I need to get value say f($X_i$, $Y_i$). How to solve this problem? I checked all the interpolation methods but nothing helped me.
Here is a numerical question format of my problem
If f(0, 1) = 3, f(2, 3) = 7, f(4, 2) = 8 and, f(2, 2) = 6, find the values at f(3, 5) using Langrange Interpolation or any other good interpolation method.
By good interpolation method I mean any method that performs interpolation with minimum errors in value approximation.
 A: For a general scattered data interpolation problem with $N$ known data points in $\mathbb{R}^d$, I recommend trying a radial basis function approach. Let the data sites be denoted $\{x_1,...,x_N\}$ with corresponding function values $\{y_j\}_{j=1}^N$. We choose a radial basis function $\varphi(r): \mathbb{R} \to \mathbb{R}$, and define $\Phi_j: \mathbb{R}^d \to \mathbb{R}$ by $\Phi_j(x) = \varphi(||x-x_j||)$. We enforce the interpolation condition by saying for each $x_k$, we want our interpolant $s(x)$ to satisfies $s(x_k) = y_k$. That is,
$$s(x_k) = \sum_{j=1}^N c_j \varphi(||x_k-x_j||)$$
This yields a system of $N$ equations in $N$ unknowns (the $c_j$ being unknown) of the form $Ac = y$, where $y$ is a vector of the known data values $y_j$, $c$ is the coefficient vector of the RBF expansion, and $A_{kj} = \varphi(||x_k-x_j||)$.
If the chosen RBF $\varphi$ is positive definite, then no matter what collection of data sites used, the matrix system $Ac=y$ will be uniquely solvable. Error approximation rates are known for certain classes of RBFs. A common one to use is a Gaussian $\varphi(r) = exp(-\epsilon r^2)$ for some $\epsilon$, a thin plate spline $\varphi(r) = r^2 \log(r)$ (conditionally positive definite, but similar idea), or a compactly supported RBF such as $\varphi(r) = (1-r)_+$.
These above methods do NOT require gridded data and they work on any scattered collection of data sites. For more information, see Holger Wendland's book "Scattered Data Approximation" or see Greg Fasshauer's website, where he has lecture notes and presentations on multivariate approximation from courses he's taught.
(See the first chapter probably)
http://www.math.iit.edu/~fass/603_handouts.html
http://www.math.iit.edu/~fass/590/notes/
Hope this helps!
A: In general, give data $(x_i,y_i)$, $i=1\ldots n$, you want $n$ functions
$f_j(x,y)$ so that the matrix $V$ with entries $v_{ij} = f_j(x_i,y_i)$ is
nonsingular.  Then to interpolate $f(x_i,y_i) = z_i$ you take
$f(x,y) = \sum_{j=1}^n c_j f_j(x,y)$ where ${\bf c} = V^{-1} {\bf z}$ (making the $c_j$ and $x_j$ into column vectors).
EDIT: In your example, $(x_i, y_i)$ are $(0,1), (2,3), (4,2), (2,2)$.  You might
try some low-degree polynomials:  $f_1(x,y) = 1$, $f_2(x,y) = x$, $f_3(x,y) = y$, $f_4(x,y) = x^2$.  Then 
$$ V = \left[ \begin {array}{cccc} f_{{1}} \left( 0,1 \right) &f_{{2}}
 \left( 0,1 \right) &f_{{3}} \left( 0,1 \right) &f_{{4}} \left( 0,1
 \right) \\ f_{{1}} \left( 2,3 \right) &f_{{2}}
 \left( 2,3 \right) &f_{{3}} \left( 2,3 \right) &f_{{4}} \left( 2,3
 \right) \\ f_{{1}} \left( 4,2 \right) &f_{{2}}
 \left( 4,2 \right) &f_{{3}} \left( 4,2 \right) &f_{{4}} \left( 4,2
 \right) \\ f_{{1}} \left( 2,2 \right) &f_{{2}}
 \left( 2,2 \right) &f_{{3}} \left( 2,2 \right) &f_{{4}} \left( 2,2
 \right) \end {array} \right] = \left[ \begin {array}{cccc} 1&0&1&0\\ 1&2&3&4
\\ 1&4&2&16\\ 1&2&2&4\end {array}
 \right] $$
which is nonsingular (its determinant is 16).
With ${\bf z} = \pmatrix{3\cr 7\cr 8\cr 6\cr}$ you get ${\bf c} = \pmatrix{2\cr 1\cr 1\cr 0\cr}$, i.e. $f(x,y) = 2 + x + y$, and in particular $f(3,5) = 10$.
