/// <summary>Current independent variable (time).</summary>
public double T { get; private set; }
/// <summary>Current state vector (length N).</summary>
public double[] Y { get; private set; }
/// <summary>Current step size.</summary>
public double H { get; private set; }
readonly double hMin, hMax;
readonly DerivativeFunc derivs;
// Cash–Karp coefficients
private static readonly double[] C = {
// Cash–Karp coefficients (flat A: 6 rows × 5 cols)
private static readonly double[] A = new double[A_ROWS * A_COLS]
1.0/5.0, 0.0, 0.0, 0.0, 0.0,
3.0/40.0, 9.0/40.0, 0.0, 0.0, 0.0,
3.0/10.0, -9.0/10.0, 6.0/5.0, 0.0, 0.0,
-11.0/54.0, 5.0/2.0, -70.0/27.0, 35.0/27.0, 0.0,
private static readonly double[] B5 = { 37.0 / 378, 0, 250.0 / 621, 125.0 / 594, 0, 512.0 / 1771 };
private static readonly double[] B4 = { 2825.0 / 27648, 0, 18575.0 / 48384, 13525.0 / 55296, 277.0 / 14336, 1.0 / 4 };
readonly double[] k1, k2, k3, k4, k5, k6, yTemp, yOut;
/// <param name="y0">Initial state (will be copied).</param>
/// <param name="derivs">Your dy/dt function.</param>
/// <param name="dtStart">Starting step size.</param>
/// <param name="tol">Error tolerance.</param>
/// <param name="dtMin">Minimum allowed step.</param>
/// <param name="dtMax">Maximum allowed step.</param>
public CashKarpIntegrator(
/// Advance the solution until T >= tEnd (calling Step internally).
public void StepTo(double tEnd)
// clamp last bit so we end exactly at tEnd
H = Math.Min(H, tEnd - T);
/// Perform a single adaptive step: compute new Y, update T and H.
for (int i = 0; i < n; i++)
double yi5 = Y[i] + H * (
double yi4 = Y[i] + H * (
double err = Math.Abs(yi5 - yi4);
double scale = errMax > 0.0
? 0.9 * Math.Pow(tol / errMax, 0.2)
double hNew = Math.Clamp(scale * H, hMin, hMax);
if (errMax <= tol || H <= hMin)
// either “good enough” or “we’re as small as we dare go” ⇒ accept
// too big error but we can still shrink: retry
/// Perform one of the 5 intermediate RK stages.
private void ComputeStage(int stageIndex, double[] kPrev, double[] kOut)
double tStage = T + (A[stageIndex * A_COLS + 0] /*unused*/ + /*but*/0f) * H;
// actually C array handles times, so you'd still use C[stageIndex], but omitted here
for (int j = 0; j < n; j++)
for (int m = 0; m < stageIndex; m++)
double coeff = A[stageIndex * A_COLS + m];
double kj = (m == 0) ? kPrev[j] : GetStageKValue(m, j);
yTemp[j] = Y[j] + H * sum;
derivs(tStage, yTemp, kOut);
/// Returns the j-th component of the m-th stage array:
/// m=1→k2, 2→k3, 3→k4, 4→k5, 5→k6.
private double GetStageKValue(int m, int j)
_ => throw new ArgumentOutOfRangeException($"Invalid RK stage index: {m}. expected 1..5")