Write it as:

$\displaystyle u(x,y)=a +bx+cy+dx^2+ey^2+fxy+gx^2y+hyx^2+ix^3+jy^3$

The following Mathematica code first generates a list {{x_1,y_1,z_1},. . .{x_n,y_n,z_n}} for the function $\displaystyle \sin(xy)$ then fits the data points to u(x,y) then superimposes the plots to obtain a quick visual indication of how well the fit was determined.

Code:

f[x_,y_]=Sin[x y];
mytable = Flatten[Table[{x, y, f[x, y]},
{x, -1, 1, 0.1}, {y, -1, 1, 0.1}],
1];
myequation = a + b*x + c*y + d*x^2 +
e*y^2 + f*x*y + g*x^2*y + h*y^2*x +
i*y^3 + j*x^3;
myconstants = FindFit[mytable,
myequation, {a, b, c, d, e, f, g, h,
i, j}, {x, y}]
pic2 = Plot3D[myequation /. myconstants,
{x, -1, 1}, {y, -1, 1},
PlotStyle -> Blue]
Show[{pic1, pic2}]