diff --git a/src/main/java/org/apache/commons/math4/ode/nonstiff/ThreeEighthesFieldIntegrator.java b/src/main/java/org/apache/commons/math4/ode/nonstiff/ThreeEighthesFieldIntegrator.java
new file mode 100644
index 000000000..c0d11b8a3
--- /dev/null
+++ b/src/main/java/org/apache/commons/math4/ode/nonstiff/ThreeEighthesFieldIntegrator.java
@@ -0,0 +1,77 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements. See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License. You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+package org.apache.commons.math4.ode.nonstiff;
+
+import org.apache.commons.math4.Field;
+import org.apache.commons.math4.RealFieldElement;
+
+/**
+ * This class implements the 3/8 fourth order Runge-Kutta
+ * integrator for Ordinary Differential Equations.
+ *
+ * This method is an explicit Runge-Kutta method, its Butcher-array
+ * is the following one :
+ *
+ * 0 | 0 0 0 0
+ * 1/3 | 1/3 0 0 0
+ * 2/3 |-1/3 1 0 0
+ * 1 | 1 -1 1 0
+ * |--------------------
+ * | 1/8 3/8 3/8 1/8
+ *
+ *
+ *
+ * @see EulerFieldIntegrator
+ * @see ClassicalRungeKuttaFieldIntegrator
+ * @see GillfieldIntegrator
+ * @see MidpointFieldIntegrator
+ * @see LutherFieldIntegrator
+ * @param the type of the field elements
+ * @since 3.6
+ */
+
+public class ThreeEighthesFieldIntegrator>
+ extends RungeKuttaFieldIntegrator {
+
+ /** Time steps Butcher array. */
+ private static final double[] STATIC_C = {
+ 1.0 / 3.0, 2.0 / 3.0, 1.0
+ };
+
+ /** Internal weights Butcher array. */
+ private static final double[][] STATIC_A = {
+ { 1.0 / 3.0 },
+ { -1.0 / 3.0, 1.0 },
+ { 1.0, -1.0, 1.0 }
+ };
+
+ /** Propagation weights Butcher array. */
+ private static final double[] STATIC_B = {
+ 1.0 / 8.0, 3.0 / 8.0, 3.0 / 8.0, 1.0 / 8.0
+ };
+
+ /** Simple constructor.
+ * Build a 3/8 integrator with the given step.
+ * @param field field to which the time and state vector elements belong
+ * @param step integration step
+ */
+ public ThreeEighthesFieldIntegrator(final Field field, final T step) {
+ super(field, "3/8", STATIC_C, STATIC_A, STATIC_B, new ThreeEighthesFieldStepInterpolator(), step);
+ }
+
+}
diff --git a/src/main/java/org/apache/commons/math4/ode/nonstiff/ThreeEighthesFieldStepInterpolator.java b/src/main/java/org/apache/commons/math4/ode/nonstiff/ThreeEighthesFieldStepInterpolator.java
new file mode 100644
index 000000000..c389f4641
--- /dev/null
+++ b/src/main/java/org/apache/commons/math4/ode/nonstiff/ThreeEighthesFieldStepInterpolator.java
@@ -0,0 +1,152 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements. See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License. You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+package org.apache.commons.math4.ode.nonstiff;
+
+import org.apache.commons.math4.RealFieldElement;
+import org.apache.commons.math4.ode.FieldEquationsMapper;
+import org.apache.commons.math4.ode.FieldODEStateAndDerivative;
+import org.apache.commons.math4.util.MathArrays;
+
+/**
+ * This class implements a step interpolator for the 3/8 fourth
+ * order Runge-Kutta integrator.
+ *
+ * This interpolator allows to compute dense output inside the last
+ * step computed. The interpolation equation is consistent with the
+ * integration scheme :
+ *
+ * - Using reference point at step start:
+ * y(tn + θ h) = y (tn)
+ * + θ (h/8) [ (8 - 15 θ + 8 θ2) y'1
+ * + 3 * (15 θ - 12 θ2) y'2
+ * + 3 θ y'3
+ * + (-3 θ + 4 θ2) y'4
+ * ]
+ *
+ * - Using reference point at step end:
+ * y(tn + θ h) = y (tn + h)
+ * - (1 - θ) (h/8) [(1 - 7 θ + 8 θ2) y'1
+ * + 3 (1 + θ - 4 θ2) y'2
+ * + 3 (1 + θ) y'3
+ * + (1 + θ + 4 θ2) y'4
+ * ]
+ *
+ *
+ *
+ *
+ * where θ belongs to [0 ; 1] and where y'1 to y'4 are the four
+ * evaluations of the derivatives already computed during the
+ * step.
+ *
+ * @see ThreeEighthesFieldIntegrator
+ * @param the type of the field elements
+ * @since 3.6
+ */
+
+class ThreeEighthesFieldStepInterpolator>
+ extends RungeKuttaFieldStepInterpolator {
+
+ /** Simple constructor.
+ * This constructor builds an instance that is not usable yet, the
+ * {@link
+ * org.apache.commons.math4.ode.sampling.AbstractFieldStepInterpolator#reinitialize}
+ * method should be called before using the instance in order to
+ * initialize the internal arrays. This constructor is used only
+ * in order to delay the initialization in some cases. The {@link
+ * RungeKuttaFieldIntegrator} class uses the prototyping design pattern
+ * to create the step interpolators by cloning an uninitialized model
+ * and later initializing the copy.
+ */
+ ThreeEighthesFieldStepInterpolator() {
+ }
+
+ /** Copy constructor.
+ * @param interpolator interpolator to copy from. The copy is a deep
+ * copy: its arrays are separated from the original arrays of the
+ * instance
+ */
+ ThreeEighthesFieldStepInterpolator(final ThreeEighthesFieldStepInterpolator interpolator) {
+ super(interpolator);
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ protected ThreeEighthesFieldStepInterpolator doCopy() {
+ return new ThreeEighthesFieldStepInterpolator(this);
+ }
+
+
+ /** {@inheritDoc} */
+ @Override
+ protected FieldODEStateAndDerivative computeInterpolatedStateAndDerivatives(final FieldEquationsMapper mapper,
+ final T time, final T theta,
+ final T oneMinusThetaH) {
+
+ final T coeffDot3 = theta.multiply(0.75);
+ final T coeffDot1 = coeffDot3.multiply(theta.multiply(4).subtract(5)).add(1);
+ final T coeffDot2 = coeffDot3.multiply(theta.multiply(-6).add(5));
+ final T coeffDot4 = coeffDot3.multiply(theta.multiply(2).subtract(1));
+ final T[] interpolatedState = MathArrays.buildArray(theta.getField(), previousState.length);
+ final T[] interpolatedDerivatives = MathArrays.buildArray(theta.getField(), previousState.length);
+
+ if ((previousState != null) && (theta.getReal() <= 0.5)) {
+ final T s = theta.multiply(h).divide(8);
+ final T fourTheta2 = theta.multiply(theta).multiply(4);
+ final T coeff1 = s.multiply(fourTheta2.multiply(2).subtract(theta.multiply(15)).add(8));
+ final T coeff2 = s.multiply(theta.multiply(5).subtract(fourTheta2)).multiply(3);
+ final T coeff3 = s.multiply(theta).multiply(3);
+ final T coeff4 = s.multiply(fourTheta2.subtract(theta.multiply(3)));
+ for (int i = 0; i < interpolatedState.length; ++i) {
+ final T yDot1 = yDotK[0][i];
+ final T yDot2 = yDotK[1][i];
+ final T yDot3 = yDotK[2][i];
+ final T yDot4 = yDotK[3][i];
+ interpolatedState[i] = previousState[i].
+ add(coeff1.multiply(yDot1)).add(coeff2.multiply(yDot2)).
+ add(coeff3.multiply(yDot3)).add(coeff4.multiply(yDot4));
+ interpolatedDerivatives[i] = coeffDot1.multiply(yDot1).add(coeffDot2.multiply(yDot2)).
+ add(coeffDot3.multiply(yDot3)).add(coeffDot4.multiply(yDot4));
+
+ }
+ } else {
+ final T s = oneMinusThetaH.divide(8);
+ final T fourTheta2 = theta.multiply(theta).multiply(4);
+ final T thetaPlus1 = theta.add(1);
+ final T coeff1 = s.multiply(fourTheta2.multiply(2).subtract(theta.multiply(7)).add(1));
+ final T coeff2 = s.multiply(thetaPlus1.subtract(fourTheta2)).multiply(3);
+ final T coeff3 = s.multiply(thetaPlus1).multiply(3);
+ final T coeff4 = s.multiply(thetaPlus1.add(fourTheta2));
+ for (int i = 0; i < interpolatedState.length; ++i) {
+ final T yDot1 = yDotK[0][i];
+ final T yDot2 = yDotK[1][i];
+ final T yDot3 = yDotK[2][i];
+ final T yDot4 = yDotK[3][i];
+ interpolatedState[i] = currentState[i].
+ subtract(coeff1.multiply(yDot1)).subtract(coeff2.multiply(yDot2)).
+ subtract(coeff3.multiply(yDot3)).subtract(coeff4.multiply(yDot4));
+ interpolatedDerivatives[i] = coeffDot1.multiply(yDot1).add(coeffDot2.multiply(yDot2)).
+ add(coeffDot3.multiply(yDot3)).add(coeffDot4.multiply(yDot4));
+
+ }
+ }
+
+ return new FieldODEStateAndDerivative(time, interpolatedState, yDotK[0]);
+
+ }
+
+}