How to optimize the fit of a response surface?

I need to find a better approximating response surface for the data acquired through a known function. Let's say I choose a function (in lack of experimental data) and fill an input/output matrix with x1,x2 (my 2 variables) and one with y (the responses). Applying

x1min=-5;
x1max=5;
x2min=-5;
x2max=5;
c2(5,1)=x1min;
c2(6,1)=x1max;
c2(1:2,1)=x1min/1.4142;
c2(3:4,1)=x1max/1.4142;
c2(7:16,1)=0.5*(x1max+x1min);
c2(7,2)=x2min;
c2(8,2)=x2max;
c2(1,2)=x2min/1.4142;
c2(3,2)=c2(1,2);
c2(2,2)=x2max/1.4142;
c2(4,2)=c2(2,2);
c2(5:6,2)=0.5*(x2max+x2min);
c2(9:16,2)=0.5*(x2max+x2min);

      for i=1:16 
          x1=c2(i,1);
          x2=c2(i,2);
          f(i)=x1.^2.*x2-2.*x2;
      end
      f=f';
      stat_out = regstats(f,c2,'quadratic');

The above script produces the "experimental data" that is the 'c2' matrix which contains the x1,x2 values in CCD Circumscribed design, the 'f' matrix containing the responses and then regstats, though regression, produces the 'b' coefficients in order to form a second-order approximating polynomial such as:

approx=b(1)+b(2).*x1(j)+b(3).*x2(i)+b(4).*x1(j).*x2(i)+b(5).*x1(j).^2+b(6).*x2(i).^2;

By plotting the original function and the approximating one I notice a huge difference between the two, so I have to find a way to optimize that fit so the response surface (resulting from 'approx') is much closer to the original function. If anyone has to offer any knowledge I would appreciate it a lot. Thanks in advance

[dpb]

The design points are wrong; I don't know what that particular design is OTOMH, but

>> c2
c2 =
 -3.5356   -3.5356
 -3.5356    3.5356
  3.5356   -3.5356
  3.5356    3.5356
 -5.0000         0
  5.0000         0
       0   -5.0000
       0    5.0000
       0         0
       0         0
       0         0
       0         0
       0         0
       0         0
       0         0
       0         0
>>

Half of your design points are the center and looking at the structure, only b(3), the coefficient of X2 is estimable; the rest are zero to machine precision.

>> format shorte
>> stat_out.beta
ans =
-4.0168e-16
 3.5527e-16
 4.2502e+00
-5.6842e-16
-2.0097e-16
 5.6523e-17
>>

[Nick Chatz]

Thnks for your reply. I understand that the problem lies in the first step of my design. What do you suggest I do to fix my points?

[Nick Chatz]

Thank you too. What should I do in terms of setting up the problem in order to get a decent aproximation?

[dpb]

OK, sorry, I spoke somewhat too quickly; the design is OK; I can't count +/- as two this morning... :) There are five levels, the problem is actually in the response data you generated creating the observed result.

You built in the correlation between x1,x2 as

y = x1^2 * x2 -2*x2 --> x2*[x1^2 - 2]

which introduced the linear correlation with x2 and makes other terms vanish.

What's the underlying problem you're trying to solve here?

[Nick Chatz]

I wish there was an underlying problem. The assignment I got is literally to approximate via RSM some "experimental data" (that occur from getting the responses of a certain function I have to choose for 2 variables, in this case

y = x1^2 * x2 -2*x2

Thus I have to produce some kind of polynomial that when plotted is an adequate approximation of the y plot.

[dpb]

OK, given that it is known the "experimental" data are as given, the problem in using standard RSM is that the model is inadequate. The quadratic RSM model has quadratic terms in Xi and interaction terms in Xij but there aren't higher-level interaction terms which is what is present here.

I had been presuming this was simply a stand-in for experimental data; not an actual evaluation of a simulation. If it is, indeed, the latter, then ground rules for design change drastically -- there is no sense in using replication points to estimate random error; there is essentially none so that's pointless.

Hence, the question then becomes one of how to distribute sample points to minimize bias presuming the model is adequate; the issue here is that a quadratic model is simply always going to be inadequate owing to not having interaction terms incorporated in it that can mimic the correlation between the coefficient of the x1 term and x2. So what real constraints are imposed? Is it mandatory to use RSM methods and not some other interpolating form? Who set such rules and why (always question authority! :) )?

ADDENDUM:

NB: Any design matrix that is symmetric around the origin is going to result in zero estimates for the terms other than X2 because the slope of the points in X1 will be zero owing to the symmetry...and that there isn't a higher-ordered term in the model to pick up the shape.