Fitting a second order polynomial with singular value decomposition

Dear all,

I'm trying to fit in IDL a function of the type:

y = a0 + a1x1 + a2x1^2 + a3x2 + a4x2^2 + a5x3 + a6x3^2

I've tried to use svdfit but I believe this procedure only allows
functions of the type y=f(x).

Any suggestion would be greatly appreciated.

Thanks and best wishes,

Joao