final T f2        = f1.multiply(eta);
            final T f3        = f2.multiply(theta);
            final T f4        = f3.multiply(eta);
            final T coeff0    = f1.multiply(a70).
                                subtract(f2.multiply(a70.subtract(1))).
                                add(f3.multiply(a70.multiply(2).subtract(1))).
                                add(f4.multiply(d0));
            final T coeff1    = time.getField().getZero();
            final T coeff2    = f1.multiply(a72).
                                subtract(f2.multiply(a72)).
                                add(f3.multiply(a72.multiply(2))).
                                add(f4.multiply(d2));
            final T coeff3    = f1.multiply(a73).
                                subtract(f2.multiply(a73)).
                                add(f3.multiply(a73.multiply(2))).
                                add(f4.multiply(d3));
            final T coeff4    = f1.multiply(a74).
                                subtract(f2.multiply(a74)).
                                add(f3.multiply(a74.multiply(2))).
                                add(f4.multiply(d4));
            final T coeff5    = f1.multiply(a75).
                                subtract(f2.multiply(a75)).
                                add(f3.multiply(a75.multiply(2))).
                                add(f4.multiply(d5));
            final T coeff6    = f4.multiply(d6).subtract(f3);
            final T coeffDot0 = a70.
                                subtract(dot2.multiply(a70.subtract(1))).
                                add(dot3.multiply(a70.multiply(2).subtract(1))).
                                add(dot4.multiply(d0));
            final T coeffDot1 = time.getField().getZero();
            final T coeffDot2 = a72.
                                subtract(dot2.multiply(a72)).
                                add(dot3.multiply(a72.multiply(2))).
                                add(dot4.multiply(d2));
            final T coeffDot3 = a73.
                                subtract(dot2.multiply(a73)).
                                add(dot3.multiply(a73.multiply(2))).
                                add(dot4.multiply(d3));
            final T coeffDot4 = a74.
                                subtract(dot2.multiply(a74)).
                                add(dot3.multiply(a74.multiply(2))).
                                add(dot4.multiply(d4));
            final T coeffDot5 = a75.
                                subtract(dot2.multiply(a75)).
                                add(dot3.multiply(a75.multiply(2))).
                                add(dot4.multiply(d5));
            final T coeffDot6 = dot4.multiply(d6).subtract(dot3);
            interpolatedState       = previousStateLinearCombination(coeff0, coeff1, coeff2, coeff3,

final T f1 = f0.multiply(eta);
            final T f2 = f1.multiply(theta);
            final T f3 = f2.multiply(eta);
            final T f4 = f3.multiply(theta);
            final T f5 = f4.multiply(eta);
            final T f6 = f5.multiply(theta);
            final T[] p = MathArrays.buildArray(time.getField(), 16);
            final T[] q = MathArrays.buildArray(time.getField(), 16);
            for (int i = 0; i < p.length; ++i) {
                p[i] =     f0.multiply(d[0][i]).
                       add(f1.multiply(d[1][i])).
                       add(f2.multiply(d[2][i])).
                       add(f3.multiply(d[3][i])).
                       add(f4.multiply(d[4][i])).
                       add(f5.multiply(d[5][i])).
                       add(f6.multiply(d[6][i]));
                q[i] =                    d[0][i].
                        add(dot1.multiply(d[1][i])).
                        add(dot2.multiply(d[2][i])).
                        add(dot3.multiply(d[3][i])).
                        add(dot4.multiply(d[4][i])).
                        add(dot5.multiply(d[5][i])).
                        add(dot6.multiply(d[6][i]));
            }
            interpolatedState       = previousStateLinearCombination(p[0], p[1], p[ 2], p[ 3], p[ 4], p[ 5], p[ 6], p[ 7],

final double f3        = f2 * theta;
            final double f4        = f3 * eta;
            final double coeff0    = f1 * A70 - f2   * (A70 - 1) + f3   * (2 * A70 - 1) + f4   * D0;
            final double coeff1    = 0;
            final double coeff2    = f1 * A72 - f2   * A72       + f3   * (2 * A72)     + f4   * D2;
            final double coeff3    = f1 * A73 - f2   * A73       + f3   * (2 * A73)     + f4   * D3;
            final double coeff4    = f1 * A74 - f2   * A74       + f3   * (2 * A74)     + f4   * D4;
            final double coeff5    = f1 * A75 - f2   * A75       + f3   * (2 * A75)     + f4   * D5;
            final double coeff6    = f4 * D6 - f3;
            final double coeffDot0 =      A70 - dot2 * (A70 - 1) + dot3 * (2 * A70 - 1) + dot4 * D0;
            final double coeffDot1 = 0;
            final double coeffDot2 =      A72 - dot2 * A72       + dot3 * (2 * A72)     + dot4 * D2;
            final double coeffDot3 =      A73 - dot2 * A73       + dot3 * (2 * A73)     + dot4 * D3;
            final double coeffDot4 =      A74 - dot2 * A74       + dot3 * (2 * A74)     + dot4 * D4;
            final double coeffDot5 =      A75 - dot2 * A75       + dot3 * (2 * A75)     + dot4 * D5;
            final double coeffDot6 = dot4 * D6 - dot3;
            interpolatedState       = previousStateLinearCombination(coeff0, coeff1, coeff2, coeff3,

final double f2 = f1 * theta;
            final double f3 = f2 * eta;
            final double f4 = f3 * theta;
            final double f5 = f4 * eta;
            final double f6 = f5 * theta;
            final double[] p = new double[16];
            final double[] q = new double[16];
            for (int i = 0; i < p.length; ++i) {
                p[i] = f0 * D[0][i] +   f1 * D[1][i] +   f2 * D[2][i] +   f3 * D[3][i] +
                       f4 * D[4][i] +   f5 * D[5][i] +   f6 * D[6][i];
                q[i] = D[0][i] + dot1 * D[1][i] + dot2 * D[2][i] + dot3 * D[3][i] +
                                 dot4 * D[4][i] + dot5 * D[5][i] + dot6 * D[6][i];
            }
            interpolatedState       = previousStateLinearCombination(p[0], p[1], p[ 2], p[ 3], p[ 4], p[ 5], p[ 6], p[ 7],

final T[][] a, final T[] c,
                                                                                final T[][] yDotK) {
        // create some internal working arrays
        final int stages = c.length + 1;
        final T[] yTmp = y0.clone();

        for (int k = 1; k < stages; ++k) {

            for (int j = 0; j < y0.length; ++j) {
                T sum = yDotK[0][j].multiply(a[k - 1][0]);
                for (int l = 1; l < k; ++l) {
                    sum = sum.add(yDotK[l][j].multiply(a[k - 1][l]));
                }
                yTmp[j] = y0[j].add(h.multiply(sum));
            }

            yDotK[k] = equations.computeDerivatives(t0.add(h.multiply(c[k - 1])), yTmp);
        }
    }

    /** Apply internal weights of Butcher array, with corresponding times. Version with real Butcher array (non-Field).
     * @param <T> the type of the field elements
     * @param equations differential equations to integrate
     * @param t0 initial time
     * @param y0 initial value of the state vector at t0
     * @param h step size
     * @param a internal weights of Butcher array
     * @param c times of Butcher array
     * @param yDotK array where to store result
     */
    static <T extends CalculusFieldElement<T>> void applyInternalButcherWeights(final FieldExpandableODE<T> equations,

}
            interpolatedState       = currentStateLinearCombination(p[0], p[1], p[ 2], p[ 3], p[ 4], p[ 5], p[ 6], p[ 7],
                                                                    p[8], p[9], p[10], p[11], p[12], p[13], p[14], p[15]);
            interpolatedDerivatives = derivativeLinearCombination(q[0], q[1], q[ 2], q[ 3], q[ 4], q[ 5], q[ 6], q[ 7],
                                                                  q[8], q[9], q[10], q[11], q[12], q[13], q[14], q[15]);
        }

        return mapper.mapStateAndDerivative(time, interpolatedState, interpolatedDerivatives);

    }

}

}
            interpolatedState       = previousStateLinearCombination(p[0], p[1], p[ 2], p[ 3], p[ 4], p[ 5], p[ 6], p[ 7],
                                                                     p[8], p[9], p[10], p[11], p[12], p[13], p[14], p[15]);
            interpolatedDerivatives = derivativeLinearCombination(q[0], q[1], q[ 2], q[ 3], q[ 4], q[ 5], q[ 6], q[ 7],
                                                                  q[8], q[9], q[10], q[11], q[12], q[13], q[14], q[15]);
        } else {
            final T f0 = oneMinusThetaH.negate();

interpolatedState       = previousStateLinearCombination(p[0], p[1], p[ 2], p[ 3], p[ 4], p[ 5], p[ 6], p[ 7],
                                                                     p[8], p[9], p[10], p[11], p[12], p[13], p[14], p[15]);
            interpolatedDerivatives = derivativeLinearCombination(q[0], q[1], q[ 2], q[ 3], q[ 4], q[ 5], q[ 6], q[ 7],
                                                                  q[8], q[9], q[10], q[11], q[12], q[13], q[14], q[15]);
        } else {

@Override
    protected FieldODEStateAndDerivative<T> initIntegration(FieldExpandableODE<T> eqn, FieldODEState<T> s0, T t) {
        if (!isUsingFieldCoefficients()) {
            realA = getRealA();
            realB = getRealB();
            realC = getRealC();
        }
        return super.initIntegration(eqn, s0, t);
    }

    /** {@inheritDoc} */
    @Override
    public FieldODEStateAndDerivative<T> integrate(final FieldExpandableODE<T> equations,
                                                   final FieldODEState<T> initialState, final T finalTime)
        throws MathIllegalArgumentException, MathIllegalStateException {

        sanityChecks(initialState, finalTime);
        setStepStart(initIntegration(equations, initialState, finalTime));
        final boolean forward = finalTime.subtract(initialState.getTime()).getReal() > 0;

        // create some internal working arrays
        final int   stages = getNumberOfStages();
        final T[][] yDotK  = MathArrays.buildArray(getField(), stages, -1);

super(field, name);
        stepsizeHelper = new StepsizeHelper(minStep, maxStep, vecAbsoluteTolerance, vecRelativeTolerance);
        resetInternalState();
    }

    /** Set the adaptive step size control parameters.
     * <p>
     * A side effect of this method is to also reset the initial
     * step so it will be automatically computed by the integrator
     * if {@link #setInitialStepSize(double) setInitialStepSize}
     * is not called by the user.
     * </p>
     * @param minimalStep minimal step (must be positive even for backward
     * integration), the last step can be smaller than this
     * @param maximalStep maximal step (must be positive even for backward
     * integration)
     * @param absoluteTolerance allowed absolute error
     * @param relativeTolerance allowed relative error
     */
    public void setStepSizeControl(final double minimalStep, final double maximalStep,
                                   final double absoluteTolerance, final double relativeTolerance) {
        stepsizeHelper = new StepsizeHelper(minimalStep, maximalStep, absoluteTolerance, relativeTolerance);
    }

    /** Set the adaptive step size control parameters.
     * <p>
     * A side effect of this method is to also reset the initial
     * step so it will be automatically computed by the integrator
     * if {@link #setInitialStepSize(double) setInitialStepSize}
     * is not called by the user.
     * </p>
     * @param minimalStep minimal step (must be positive even for backward
     * integration), the last step can be smaller than this
     * @param maximalStep maximal step (must be positive even for backward
     * integration)
     * @param absoluteTolerance allowed absolute error
     * @param relativeTolerance allowed relative error
     */
    public void setStepSizeControl(final double minimalStep, final double maximalStep,
                                   final double[] absoluteTolerance,
                                   final double[] relativeTolerance) {
        stepsizeHelper = new StepsizeHelper(minimalStep, maximalStep, absoluteTolerance, relativeTolerance);
    }

    /** Get the stepsize helper.
     * @return stepsize helper
     * @since 2.0
     */
    protected StepsizeHelper getStepSizeHelper() {
        return stepsizeHelper;
    }

    /** Set the initial step size.
     * <p>This method allows the user to specify an initial positive
     * step size instead of letting the integrator guess it by
     * itself. If this method is not called before integration is
     * started, the initial step size will be estimated by the
     * integrator.</p>
     * @param initialStepSize initial step size to use (must be positive even
     * for backward integration ; providing a negative value or a value
     * outside of the min/max step interval will lead the integrator to
     * ignore the value and compute the initial step size by itself)
     */
    public void setInitialStepSize(final double initialStepSize) {
        stepsizeHelper.setInitialStepSize(initialStepSize);
    }

    /** {@inheritDoc} */
    @Override
    protected void sanityChecks(final FieldODEState<T> initialState, final T t)

final T h, final T[] b) {
        final T[] y = y0.clone();
        final int stages = b.length;
        for (int j = 0; j < y0.length; ++j) {
            T sum = yDotK[0][j].multiply(b[0]);
            for (int l = 1; l < stages; ++l) {
                sum = sum.add(yDotK[l][j].multiply(b[l]));
            }
            y[j] = y[j].add(h.multiply(sum));
        }
        return y;
    }

EquationsMapper(final EquationsMapper mapper, final int dimension) {
        final int index = (mapper == null) ? 0 : mapper.getNumberOfEquations();
        this.start = new int[index + 2];
        if (mapper == null) {
            start[0] = 0;
        } else {
            System.arraycopy(mapper.start, 0, start, 0, index + 1);
        }
        start[index + 1] = start[index] + dimension;
    }

    /** Get the number of equations mapped.
     * @return number of equations mapped
     */
    public int getNumberOfEquations() {
        return start.length - 1;
    }

    /** Return the dimension of the complete set of equations.
     * <p>
     * The complete set of equations correspond to the primary set plus all secondary sets.
     * </p>
     * @return dimension of the complete set of equations
     */
    public int getTotalDimension() {
        return start[start.length - 1];
    }

    /** Map flat arrays to a state and derivative.
     * @param t time
     * @param y state array to map, including primary and secondary components
     * @param yDot state derivative array to map, including primary and secondary components
     * @return mapped state
     * @exception MathIllegalArgumentException if an array does not match total dimension
     */
    public ODEStateAndDerivative mapStateAndDerivative(final double t, final double[] y, final double[] yDot)

} else if (g < 0) {
                            // initialize as if previous root (i.e. forward one) was an ignored increasing event
                            return Transformer.MIN;
                        } else {
                            // we are exactly at a root, we don't know if it is an increasing
                            // or a decreasing event, we remain in uninitialized state
                            return Transformer.UNINITIALIZED;
                        }
                    case PLUS  :
                        if (g <= 0) {
                            // we have crossed the zero line on an ignored increasing event,
                            // we must change the transformer
                            return Transformer.MAX;
                        } else {
                            // we are still in the same status
                            return previous;
                        }
                    case MINUS :
                        if (g <= 0) {
                            // we have crossed the zero line on an ignored increasing event,
                            // we must change the transformer
                            return Transformer.MIN;
                        } else {
                            // we are still in the same status
                            return previous;
                        }
                    case MIN   :
                        if (g >= 0) {
                            // we have crossed the zero line on a triggered decreasing event,
                            // we must change the transformer
                            return Transformer.PLUS;
                        } else {
                            // we are still in the same status
                            return previous;
                        }
                    case MAX   :
                        if (g >= 0) {
                            // we have crossed the zero line on a triggered decreasing event,
                            // we must change the transformer
                            return Transformer.MINUS;
                        } else {
                            // we are still in the same status
                            return previous;
                        }
                    default    :
                        // this should never happen
                        throw MathRuntimeException.createInternalError();
                }
            }

public void insertEquationData(final int index, double[] equationData, double[] complete)
        throws MathIllegalArgumentException {
        checkIndex(index);
        final int begin     = start[index];
        final int end       = start[index + 1];
        final int dimension = end - begin;
        if (complete.length < end) {
            throw new MathIllegalArgumentException(LocalizedCoreFormats.DIMENSIONS_MISMATCH,
                                                   complete.length, end);
        }
        if (equationData.length != dimension) {
            throw new MathIllegalArgumentException(LocalizedCoreFormats.DIMENSIONS_MISMATCH,
                                                   equationData.length, dimension);
        }
        System.arraycopy(equationData, 0, complete, begin, dimension);
    }

    /** Check equation index.
     * @param index index of the equation, must be between 0 included and
     * {@link #getNumberOfEquations()} (excluded)
     * @exception MathIllegalArgumentException if index is out of range
     */
    private void checkIndex(final int index) throws MathIllegalArgumentException {

c[3] = FieldExplicitRungeKuttaIntegrator.fraction(getField(),8,  9);
        c[4] = getField().getOne();
        c[5] = getField().getOne();
        return c;
    }

    /** {@inheritDoc} */
    @Override
    public T[][] getA() {
        final T[][] a = MathArrays.buildArray(getField(), 6, -1);
        for (int i = 0; i < a.length; ++i) {
            a[i] = MathArrays.buildArray(getField(), i + 1);
        }
        a[0][0] = FieldExplicitRungeKuttaIntegrator.fraction(getField(),      1,     5);

yDDotOnScale += ratioDotDot * ratioDotDot;
        }
        yDDotOnScale = FastMath.sqrt(yDDotOnScale) / h;

        // step size is computed such that
        // h^order * max (||y'/tol||, ||y''/tol||) = 0.01
        final double maxInv2 = FastMath.max(FastMath.sqrt(yDotOnScale2), yDDotOnScale);
        final double h1 = (maxInv2 < 1.0e-15) ?
                          FastMath.max(1.0e-6, 0.001 * FastMath.abs(h)) :
                          FastMath.pow(0.01 / maxInv2, 1.0 / order);
        h = FastMath.min(100.0 * FastMath.abs(h), h1);
        h = FastMath.max(h, 1.0e-12 * FastMath.abs(state0.getTime().getReal()));  // avoids cancellation when computing t1 - t0

if (isUsingFieldCoefficients()) {
                    FieldExplicitRungeKuttaIntegrator.applyInternalButcherWeights(getEquations(),
                            getStepStart().getTime(), y, getStepSize(), a, c, yDotK);
                    yTmp = FieldExplicitRungeKuttaIntegrator.applyExternalButcherWeights(y, yDotK, getStepSize(), b);
                } else {
                    FieldExplicitRungeKuttaIntegrator.applyInternalButcherWeights(getEquations(),
                            getStepStart().getTime(), y, getStepSize(), realA, realC, yDotK);
                    yTmp = FieldExplicitRungeKuttaIntegrator.applyExternalButcherWeights(y, yDotK, getStepSize(), realB);
                }

                incrementEvaluations(stages - 1);

final Transformer next     = filter.selectTransformer(previous, rawG, forward);
                if (next != previous) {
                    // there is a root somewhere between extremeT and t.
                    // the new transformer is valid for t (this is how we have just computed
                    // it above), but it is in fact valid on both sides of the root, so
                    // it was already valid before t and even up to previous time. We store
                    // the switch at extremeT for safety, to ensure the previous transformer
                    // is not applied too close of the root
                    System.arraycopy(updates,      0, updates,      1, updates.length - 1);
                    System.arraycopy(transformers, 0, transformers, 1, transformers.length - 1);
                    updates[0]      = extremeT;
                    transformers[0] = next;
                }

                extremeT = state.getTime();

                // apply the transform
                return next.transformed(rawG);

            } else {
                // we are in the middle of the history

                // select the transformer
                for (int i = 0; i < updates.length - 1; ++i) {
                    if (state.getTime() <= updates[i]) {