polynomial fitting

I'm trying to find a way to fit a polynomial of degree n in x  and y to
a set of x, y, and z data that I have and obtain the coefficients for
the terms of the fitted polynomial.  However, when I try to use the
surf.ls function I'm getting odd results.

x <- seq(0, 10, length=50)
y <- x
f <- function (x, y) {x^2 + y}  
library(spatial)
test <- data.frame(x=x, y=y, z=f(x, y))
test.kr <- surf.ls(2, test)
test.kr$beta

[1] 0 0 0 0 0 0

When I try the example from the help I get:

     library(MASS)
     data(topo, package="MASS")
     topo.kr <- surf.ls(2, topo)
topo.kr$beta

[1] 801.217617 -11.018887  68.229148 -73.992968   3.343573   8.342717

Why is my test data causing problems?  Also it seems that the beta
attribute from the object returned the surf.ls correspond with the terms
of the fitted polynomial.  If this is correct, in what order are the
coefficients for the fitted polynomial given?  Finally, the R source
code for the spatial library indicate that both surf.ls and surf.gls are
limited to polynomials of degree 6 or below.  Is there another function
that will work with higher order polynomials.  

I'm working with R 1.7.0 using the binary rpm for Redhat 7.3 from CRAN.

version

_                
platform i686-pc-linux-gnu
arch     i686             
os       linux-gnu        
system   i686, linux-gnu  
status                    
major    1                
minor    7.0              
year     2003             
month    04               
day      16               
language R