public static void polcof(final double[] xa, final double[] ya, final double[] cof) 
            throws NRException {
        // Given arrays xa[0..n-1] and ya[0..n-1] containing a tabulated
        // function yai D f.xai /, this routine returns an array of
        // coefficients cof[0..n-1], such that yai D
        // Pn1
        // jD0 cofj xa
        // j
        // i .

        int k, j, i, n = xa.length;
        double xmin;
        final double[] x = doub_vec(n), y = doub_vec(n);
        for (j = 0; j < n; j++) {
            x[j] = xa[j];
            y[j] = ya[j];
        }
        for (j = 0; j < n; j++) { // Fill a temporary vector whose size
            // decreases as each coefficient is
            // found.
            final double[] x_t = doub_vec(n - j), y_t = doub_vec(n - j);
            for (k = 0; k < n - j; k++) {
                x_t[k] = x[k];
                y_t[k] = y[k];
            }
            Poly_interp interp = new Poly_interp(x, y, n - j);
            cof[j] = interp.rawinterp(0, 0.); // Extrapolate to x D 0.
            xmin = 1.0e99;
            k = -1;
            for (i = 0; i < n - j; i++) { // Find the remaining xi of smallest
                if (abs(x[i]) < xmin) { // absolute value
                    xmin = abs(x[i]);
                    k = i;
                }
                if (x[i] != 0.0) // (meanwhile reducing all the terms)
                    y[i] = (y[i] - cof[j]) / x[i];
            }
            for (i = k + 1; i < n - j; i++) { // and eliminate it.
                y[i - 1] = y[i];
                x[i - 1] = x[i];
            }
        }
    }