How to approximate unknown two-dimensional function? I have a surface defined by values on a two-dimensional grid. I would like to find an approximate function which would give me a value for any arbitrary point within certain range of xs and ys.
My general idea is to construct some sort of polynomial, and then tweak it's coefficients by using some sort of evolutionary algorithm until the polynomial behaves as I want it to. But how should my polynomial look in general?

The optimal method depends on the structure of your grid and your problem in general. 
The first question you should ask yourself: Does your problem require that the resulting function exactly coincides with the given values at the grid points? This approach would be called interpolation, e.g. polynomial interpolation as suggested by the previous answer. Downsides of that approach are the following:
1) the resulting function is going to have strong oscillations between the gridpoints if the gridpoints lie close to each other.
2) If you are interested in the values outside the grid, you won't be able to use the interpolating polynomial at all.
If these are serious restrictions to your use case, you should consider another approach where you approximate instead of interpolate (usually using a least-squares or $L^2$ approach). The most prominent approach here is (linear) regression.
In the end, what method works better depends on where your data is coming from and what you are trying to infer from it. In any case, for both of these approaches there are many free software implementations available. 
