Suppose we have a one-variable function f(x) whose values are known at points which may be unevenly spaced:
 

We can do an exact fit using the method of Undetermined Coefficients. If we approximate the function by a linear combination of the first N_x+1 members of a certain complete set of functions {G_i}, using a series expansion of the form:

N_x
ε A_i G_i(x)   ;   G₀(x)=1
i=0

where the A_i's are constants and N_x+1 is the number of terms in x.  For an exact fit, N_x=P_x, otherwise N_x<P_x.

If we use x_C as center of expansion, the approximating function's value at the m_th grid point is:

N_x
ε A_i G_i(X_m-x_C)
i=0

To do an exact fit, at each grid point we equate the value of the function at that grid point to the value of the approximation at the same grid point:

               N_x
f(x_m) = ε A_i G_i(X_m-x_C)
       i=0

This way we set up a system of linear equations from which the coefficients A_i can be easily found.

If we have a 2-variable tabulated function f(x,y) whose values are known at the grid points of a rectangular grid which may be unevenly spaced, such as

which has P_x intervals in the x-direction and P_y intervals in the y-direction, we can approximate the function by a linear combination of the first N_x+1 and N_y+1 members of certain complete sets of functions {G_i} and {H_j} using a series expansion of the form:

N_x        N_y
ε    ε   A_ij G_i(x) H_j(y)  ;  G₀(x)=1  ;  H₀(y)=1
i=0   j=0

which has N_x+1 terms in x and N_y+1 terms in y.  For an exact fit, N_x=P_x and N_y=P_y, otherwise N_x<P_x and/or N_y<P_y.

If we use (x_C,y_C) as center of expansion, the approximating function's value at the (x_m,y_n) grid point is:

N_x        N_y
ε    ε   A_ij G_i(x_m-x_C) H_j(y_n-y_C)
i=0   j=0

To do an exact fit, at each grid point we equate the value of the function at that grid point to the value of the approximation at the same grid point:

                     N_x       N_y
f(x_m,y_n) = ε    ε   A_ij G_i(x_m-x_C) H_j(y_n-y_C)
          i=0   j=0

This way we set up a system of linear equations from which the coefficients A_ij can be easily found.

Extension to more than 2 variables is obvious.

The sets of functions used here, that is, {G_i} and {H_j}, are:

1, (x-x_C), (x-x_C)², (x-x_C)³, (x-x_C)⁴,...

and

        2p(x-x_C)         2p(x-x_C)         2*2p(x-x_C)         2*2p(x-x_C)
1, SIN( ———————— ), COS( ————————— ), SIN( —————————— ), COS( ——————————— ),...
          PER              PER                PER                PER

where x_C is the center of expansion and PER is the period.

The use of finite power and trigonometric series is partly justified by Weierstrass' Approximation Theorems:

1) Any function which is continuous in an interval may be approximated in that interval, to any degree of accuracy, by a finite power series.

2) Any continuous periodic function may be approximated, to any degree of accuracy, by a finite trigonometric series.

Naturally, if the function itself is not a finite power or trigonometric series, then we need more and more terms for more and more accuracy.  A power or trigonometric series with infinite number of terms may yield the exact value of the function, as is the case with Taylor and Fourier series.  Developing a Fourier or Taylor series, however, may not be very easy.  A Fourier series requires the evaluation of the Fourier Coefficients, which are integrals which in turn may need numerical evaluation.  A Taylor expansion requires you to find all the function's derivatives and this may be a difficult task.  (Curiously, continuous nondifferentiable functions exist, that is, there are continuous curves with no tangent at any point.)

Thus, the Weierstrass' Theorems make this approach acceptable if the solution is known to be continuous in the interval of interest.  If necessary, the problem should be solved piecewise, each time choosing an interval in which the solution and its derivatives are known to be continuous and sufficiently smooth (with a limited number of “wiggles”).