public void computeCoefficients(final int mu, final double h) {

    if ((polynomials == null) || (polynomials.length <= (mu + 4))) {
      resetTables(mu + 4);
    }

    currentDegree = mu + 4;

    for (int i = 0; i < currentState.length; ++i) {

      final double yp0   = h * y0Dot[i];
      final double yp1   = h * y1Dot[i];
      final double ydiff = y1[i] - currentState[i];
      final double aspl  = ydiff - yp1;
      final double bspl  = yp0 - ydiff;

      polynomials[0][i] = currentState[i];
      polynomials[1][i] = ydiff;
      polynomials[2][i] = aspl;
      polynomials[3][i] = bspl;

      if (mu < 0) {
        return;
      }

      // compute the remaining coefficients
      final double ph0 = 0.5 * (currentState[i] + y1[i]) + 0.125 * (aspl + bspl);
      polynomials[4][i] = 16 * (yMidDots[0][i] - ph0);

      if (mu > 0) {
        final double ph1 = ydiff + 0.25 * (aspl - bspl);
        polynomials[5][i] = 16 * (yMidDots[1][i] - ph1);

        if (mu > 1) {
          final double ph2 = yp1 - yp0;
          polynomials[6][i] = 16 * (yMidDots[2][i] - ph2 + polynomials[4][i]);

          if (mu > 2) {
            final double ph3 = 6 * (bspl - aspl);
            polynomials[7][i] = 16 * (yMidDots[3][i] - ph3 + 3 * polynomials[5][i]);

            for (int j = 4; j <= mu; ++j) {
              final double fac1 = 0.5 * j * (j - 1);
              final double fac2 = 2 * fac1 * (j - 2) * (j - 3);
              polynomials[j+4][i] =
                  16 * (yMidDots[j][i] + fac1 * polynomials[j+2][i] - fac2 * polynomials[j][i]);
            }

          }
        }
      }
    }

  }