Quadratic surface fitting

At one time I came across the need to do something like, say I had these 9 z values like heights on a map...

For points (-1,-1),...(1,1)

What I want is some function in x and y that matches these heights. Today I thought of the way to do that. First, interpolate and find the three quadratic equations one for each row. Maple can do this with the polynomial interpolation command.
The 3 equations end up being:

Now the trick I use to link these 3 equations together smoothly is I think, I would like the first equation to be valid when y=-1, and not when y=0 or y=1. I can interpolate the 3 points (-1,1),(0,0),(1,0) and multiply the first equation by that curve, and similarly for the other two:

Then I can make the whole equation by adding the 3 equations in x each multiplied by their corresponding y equation:

it simplifies to:

(You might have to click and zoom in to have a look)

The graph is:

Examining it is does go through the 9 points smoothly. I think a similar method will work for any number of points producing a corresponding equation in x and y, the highest power being the number of points in one row minus one.