Wrapper function to do polynomial fits with gsl

I have two arrays of double, say x and y containing some physical data.
I would like to find the best-fitting polynomial function at a given order:

\$y=pol_n(x)=c_0 + c_1x + c_2x^2 + ... + c_n * x^n\$

and return the coefficients (\$c\$) in an std::vector. There are no constraints on the function interface. The GNU Scientific Library, gsl, contains some handy functions to properly do the maths, however those functions are in plain old C. I have encapsulated them in this piece of C++ code:

std::vector gsl_polynomial_fit(const double * const data_x,
                                       const double * const data_y,
                                       const int n, const int order,
                                       double & chisq) {
  gsl_vector *y, *c;
  gsl_matrix *X, *cov;
  y = gsl_vector_alloc (n);
  c = gsl_vector_alloc (order+1);
  X   = gsl_matrix_alloc (n, order+1);
  cov = gsl_matrix_alloc (order+1, order+1);

  for (int i = 0; i  vc;
  for (int i = 0; i < order+1; i++) {
    vc.push_back(gsl_vector_get(c,i));
  }

  gsl_vector_free (y);
  gsl_vector_free (c);
  gsl_matrix_free (X);
  gsl_matrix_free (cov);

  return vc;
}

Although it is pretty much self-contained it doesn't look very safe. With speed and code clarity in mind, what could be done to improve it?

One issue is numerical stability. Although the GSL uses stable routines, the problem of fitting a polynomial in the monomial basis is so ill-conditioned (especially for large ranges in x and large (eg 50) degree), that you risk getting significant errors. Unless you really, really need the coefficients in this basis then I think you should use another basis, for example Legendre or Chebycheff polynomials. This will mean that any code that uses the coefficients will need to know how to do so, but the evaluation is in fact almost as efficient as evaluation in the polynomial basis.

StackExchange Code Review Q#71300, answer score: 2