SciML · ChrisRackauckas · Aug 28, 2026 · Aug 28, 2026 · Aug 28, 2026 · Aug 28, 2026
diff --git a/markdown/DAE/ChemicalAkzoNobel.md b/markdown/DAE/ChemicalAkzoNobel.md
diff --git a/markdown/DAE/LinearDAE.md b/markdown/DAE/LinearDAE.md
diff --git a/markdown/DAE/NANDGateProblem.md b/markdown/DAE/NANDGateProblem.md
diff --git a/markdown/DAE/OregoDAE.md b/markdown/DAE/OregoDAE.md
diff --git a/markdown/DAE/ROBERDAE.md b/markdown/DAE/ROBERDAE.md
diff --git a/markdown/DAE/TransistorAmplifier.md b/markdown/DAE/TransistorAmplifier.md
diff --git a/markdown/DAE/andrews_mechanism.md b/markdown/DAE/andrews_mechanism.md
diff --git a/markdown/DAE/caraxis.md b/markdown/DAE/caraxis.md
diff --git a/markdown/DAE/charge_pump.md b/markdown/DAE/charge_pump.md
diff --git a/markdown/DAE/fekete.md b/markdown/DAE/fekete.md
diff --git a/markdown/DAE/figures/ChemicalAkzoNobel_10_1.png b/markdown/DAE/figures/ChemicalAkzoNobel_10_1.png
diff --git a/markdown/DAE/figures/ChemicalAkzoNobel_2_1.png b/markdown/DAE/figures/ChemicalAkzoNobel_2_1.png
diff --git a/markdown/DAE/figures/ChemicalAkzoNobel_3_1.png b/markdown/DAE/figures/ChemicalAkzoNobel_3_1.png
diff --git a/markdown/DAE/figures/ChemicalAkzoNobel_4_1.png b/markdown/DAE/figures/ChemicalAkzoNobel_4_1.png
diff --git a/markdown/DAE/figures/ChemicalAkzoNobel_5_1.png b/markdown/DAE/figures/ChemicalAkzoNobel_5_1.png
diff --git a/markdown/DAE/figures/ChemicalAkzoNobel_6_1.png b/markdown/DAE/figures/ChemicalAkzoNobel_6_1.png
diff --git a/markdown/DAE/figures/ChemicalAkzoNobel_7_1.png b/markdown/DAE/figures/ChemicalAkzoNobel_7_1.png
diff --git a/markdown/DAE/figures/ChemicalAkzoNobel_8_1.png b/markdown/DAE/figures/ChemicalAkzoNobel_8_1.png
diff --git a/markdown/DAE/figures/ChemicalAkzoNobel_9_1.png b/markdown/DAE/figures/ChemicalAkzoNobel_9_1.png
diff --git a/markdown/DAE/figures/LinearDAE_10_1.png b/markdown/DAE/figures/LinearDAE_10_1.png
diff --git a/markdown/DAE/figures/LinearDAE_11_1.png b/markdown/DAE/figures/LinearDAE_11_1.png
diff --git a/markdown/DAE/figures/LinearDAE_12_1.png b/markdown/DAE/figures/LinearDAE_12_1.png
diff --git a/markdown/DAE/figures/LinearDAE_14_1.png b/markdown/DAE/figures/LinearDAE_14_1.png
diff --git a/markdown/DAE/figures/LinearDAE_15_1.png b/markdown/DAE/figures/LinearDAE_15_1.png
diff --git a/markdown/DAE/figures/LinearDAE_16_1.png b/markdown/DAE/figures/LinearDAE_16_1.png
diff --git a/markdown/DAE/figures/LinearDAE_17_1.png b/markdown/DAE/figures/LinearDAE_17_1.png
diff --git a/markdown/DAE/figures/LinearDAE_17_2.png b/markdown/DAE/figures/LinearDAE_17_2.png
diff --git a/markdown/DAE/figures/LinearDAE_17_3.png b/markdown/DAE/figures/LinearDAE_17_3.png
diff --git a/markdown/DAE/figures/LinearDAE_17_4.png b/markdown/DAE/figures/LinearDAE_17_4.png
diff --git a/markdown/DAE/figures/LinearDAE_17_5.png b/markdown/DAE/figures/LinearDAE_17_5.png
diff --git a/markdown/DAE/figures/LinearDAE_17_6.png b/markdown/DAE/figures/LinearDAE_17_6.png
diff --git a/markdown/DAE/figures/LinearDAE_18_1.png b/markdown/DAE/figures/LinearDAE_18_1.png
diff --git a/markdown/DAE/figures/LinearDAE_18_2.png b/markdown/DAE/figures/LinearDAE_18_2.png
diff --git a/markdown/DAE/figures/LinearDAE_18_3.png b/markdown/DAE/figures/LinearDAE_18_3.png
diff --git a/markdown/DAE/figures/LinearDAE_18_4.png b/markdown/DAE/figures/LinearDAE_18_4.png
diff --git a/markdown/DAE/figures/LinearDAE_18_5.png b/markdown/DAE/figures/LinearDAE_18_5.png
diff --git a/markdown/DAE/figures/LinearDAE_18_6.png b/markdown/DAE/figures/LinearDAE_18_6.png
diff --git a/markdown/DAE/figures/LinearDAE_19_1.png b/markdown/DAE/figures/LinearDAE_19_1.png
diff --git a/markdown/DAE/figures/LinearDAE_20_1.png b/markdown/DAE/figures/LinearDAE_20_1.png
diff --git a/markdown/DAE/figures/LinearDAE_20_2.png b/markdown/DAE/figures/LinearDAE_20_2.png
diff --git a/markdown/DAE/figures/LinearDAE_20_3.png b/markdown/DAE/figures/LinearDAE_20_3.png
diff --git a/markdown/DAE/figures/LinearDAE_20_4.png b/markdown/DAE/figures/LinearDAE_20_4.png
diff --git a/markdown/DAE/figures/LinearDAE_20_5.png b/markdown/DAE/figures/LinearDAE_20_5.png
diff --git a/markdown/DAE/figures/LinearDAE_20_6.png b/markdown/DAE/figures/LinearDAE_20_6.png
diff --git a/markdown/DAE/figures/LinearDAE_9_1.png b/markdown/DAE/figures/LinearDAE_9_1.png
diff --git a/markdown/DAE/figures/NANDGateProblem_10_1.png b/markdown/DAE/figures/NANDGateProblem_10_1.png
diff --git a/markdown/DAE/figures/NANDGateProblem_11_1.png b/markdown/DAE/figures/NANDGateProblem_11_1.png
diff --git a/markdown/DAE/figures/NANDGateProblem_12_1.png b/markdown/DAE/figures/NANDGateProblem_12_1.png
diff --git a/markdown/DAE/figures/NANDGateProblem_13_1.png b/markdown/DAE/figures/NANDGateProblem_13_1.png
diff --git a/markdown/DAE/figures/NANDGateProblem_14_1.png b/markdown/DAE/figures/NANDGateProblem_14_1.png
diff --git a/markdown/DAE/figures/NANDGateProblem_8_1.png b/markdown/DAE/figures/NANDGateProblem_8_1.png
diff --git a/markdown/DAE/figures/NANDGateProblem_9_1.png b/markdown/DAE/figures/NANDGateProblem_9_1.png
diff --git a/markdown/DAE/figures/OregoDAE_10_1.png b/markdown/DAE/figures/OregoDAE_10_1.png
diff --git a/markdown/DAE/figures/OregoDAE_2_1.png b/markdown/DAE/figures/OregoDAE_2_1.png
diff --git a/markdown/DAE/figures/OregoDAE_3_1.png b/markdown/DAE/figures/OregoDAE_3_1.png
diff --git a/markdown/DAE/figures/OregoDAE_4_1.png b/markdown/DAE/figures/OregoDAE_4_1.png
diff --git a/markdown/DAE/figures/OregoDAE_5_1.png b/markdown/DAE/figures/OregoDAE_5_1.png
diff --git a/markdown/DAE/figures/OregoDAE_6_1.png b/markdown/DAE/figures/OregoDAE_6_1.png
diff --git a/markdown/DAE/figures/OregoDAE_7_1.png b/markdown/DAE/figures/OregoDAE_7_1.png
diff --git a/markdown/DAE/figures/OregoDAE_8_1.png b/markdown/DAE/figures/OregoDAE_8_1.png
diff --git a/markdown/DAE/figures/OregoDAE_9_1.png b/markdown/DAE/figures/OregoDAE_9_1.png
diff --git a/markdown/DAE/figures/ROBERDAE_10_1.png b/markdown/DAE/figures/ROBERDAE_10_1.png
diff --git a/markdown/DAE/figures/ROBERDAE_2_1.png b/markdown/DAE/figures/ROBERDAE_2_1.png
diff --git a/markdown/DAE/figures/ROBERDAE_3_1.png b/markdown/DAE/figures/ROBERDAE_3_1.png
diff --git a/markdown/DAE/figures/ROBERDAE_4_1.png b/markdown/DAE/figures/ROBERDAE_4_1.png
diff --git a/markdown/DAE/figures/ROBERDAE_5_1.png b/markdown/DAE/figures/ROBERDAE_5_1.png
diff --git a/markdown/DAE/figures/ROBERDAE_6_1.png b/markdown/DAE/figures/ROBERDAE_6_1.png
diff --git a/markdown/DAE/figures/ROBERDAE_7_1.png b/markdown/DAE/figures/ROBERDAE_7_1.png
diff --git a/markdown/DAE/figures/ROBERDAE_8_1.png b/markdown/DAE/figures/ROBERDAE_8_1.png
diff --git a/markdown/DAE/figures/ROBERDAE_9_1.png b/markdown/DAE/figures/ROBERDAE_9_1.png
diff --git a/markdown/DAE/figures/TransistorAmplifier_2_1.png b/markdown/DAE/figures/TransistorAmplifier_2_1.png
diff --git a/markdown/DAE/figures/TransistorAmplifier_3_1.png b/markdown/DAE/figures/TransistorAmplifier_3_1.png
diff --git a/markdown/DAE/figures/TransistorAmplifier_4_1.png b/markdown/DAE/figures/TransistorAmplifier_4_1.png
diff --git a/markdown/DAE/figures/TransistorAmplifier_5_1.png b/markdown/DAE/figures/TransistorAmplifier_5_1.png
diff --git a/markdown/DAE/figures/TransistorAmplifier_6_1.png b/markdown/DAE/figures/TransistorAmplifier_6_1.png
diff --git a/markdown/DAE/figures/TransistorAmplifier_7_1.png b/markdown/DAE/figures/TransistorAmplifier_7_1.png
diff --git a/markdown/DAE/figures/TransistorAmplifier_8_1.png b/markdown/DAE/figures/TransistorAmplifier_8_1.png
diff --git a/markdown/DAE/figures/TransistorAmplifier_9_1.png b/markdown/DAE/figures/TransistorAmplifier_9_1.png
diff --git a/markdown/DAE/figures/andrews_mechanism_20_1.png b/markdown/DAE/figures/andrews_mechanism_20_1.png
diff --git a/markdown/DAE/figures/blank b/markdown/DAE/figures/blank
diff --git a/markdown/DAE/figures/caraxis_10_1.png b/markdown/DAE/figures/caraxis_10_1.png
diff --git a/markdown/DAE/figures/caraxis_11_1.png b/markdown/DAE/figures/caraxis_11_1.png
diff --git a/markdown/DAE/figures/charge_pump_15_1.png b/markdown/DAE/figures/charge_pump_15_1.png
diff --git a/markdown/DAE/figures/charge_pump_16_1.png b/markdown/DAE/figures/charge_pump_16_1.png
diff --git a/markdown/DAE/figures/fekete_13_1.png b/markdown/DAE/figures/fekete_13_1.png
diff --git a/markdown/DAE/figures/fekete_15_1.png b/markdown/DAE/figures/fekete_15_1.png
diff --git a/markdown/DAE/figures/fekete_16_1.png b/markdown/DAE/figures/fekete_16_1.png
diff --git a/markdown/DAE/figures/fekete_17_1.png b/markdown/DAE/figures/fekete_17_1.png
diff --git a/markdown/DAE/figures/fekete_18_1.png b/markdown/DAE/figures/fekete_18_1.png
diff --git a/markdown/DAE/figures/fekete_19_1.png b/markdown/DAE/figures/fekete_19_1.png
diff --git a/markdown/DAE/figures/fekete_20_1.png b/markdown/DAE/figures/fekete_20_1.png
diff --git a/markdown/DAE/figures/slider_crank_10_1.png b/markdown/DAE/figures/slider_crank_10_1.png
diff --git a/markdown/DAE/figures/slider_crank_13_1.png b/markdown/DAE/figures/slider_crank_13_1.png
diff --git a/markdown/DAE/figures/slider_crank_14_1.png b/markdown/DAE/figures/slider_crank_14_1.png
diff --git a/markdown/DAE/figures/slider_crank_15_1.png b/markdown/DAE/figures/slider_crank_15_1.png
diff --git a/markdown/DAE/figures/two_bit_adder_14_1.png b/markdown/DAE/figures/two_bit_adder_14_1.png
diff --git a/markdown/DAE/figures/two_bit_adder_16_1.png b/markdown/DAE/figures/two_bit_adder_16_1.png
diff --git a/markdown/DAE/figures/two_bit_adder_17_1.png b/markdown/DAE/figures/two_bit_adder_17_1.png
diff --git a/markdown/DAE/figures/two_bit_adder_20_1.png b/markdown/DAE/figures/two_bit_adder_20_1.png
diff --git a/markdown/DAE/figures/two_bit_adder_21_1.png b/markdown/DAE/figures/two_bit_adder_21_1.png
diff --git a/markdown/DAE/figures/two_bit_adder_22_1.png b/markdown/DAE/figures/two_bit_adder_22_1.png
diff --git a/markdown/DAE/figures/water_tube_12_1.png b/markdown/DAE/figures/water_tube_12_1.png
diff --git a/markdown/DAE/figures/water_tube_13_1.png b/markdown/DAE/figures/water_tube_13_1.png
diff --git a/markdown/DAE/figures/water_tube_14_1.png b/markdown/DAE/figures/water_tube_14_1.png
diff --git a/markdown/DAE/figures/water_tube_15_1.png b/markdown/DAE/figures/water_tube_15_1.png
diff --git a/markdown/DAE/figures/water_tube_16_1.png b/markdown/DAE/figures/water_tube_16_1.png
diff --git a/markdown/DAE/figures/water_tube_17_1.png b/markdown/DAE/figures/water_tube_17_1.png
diff --git a/markdown/DAE/figures/water_tube_18_1.png b/markdown/DAE/figures/water_tube_18_1.png
diff --git a/markdown/DAE/figures/water_tube_20_1.png b/markdown/DAE/figures/water_tube_20_1.png
diff --git a/markdown/DAE/figures/water_tube_21_1.png b/markdown/DAE/figures/water_tube_21_1.png
diff --git a/markdown/DAE/figures/water_tube_22_1.png b/markdown/DAE/figures/water_tube_22_1.png
diff --git a/markdown/DAE/figures/water_tube_23_1.png b/markdown/DAE/figures/water_tube_23_1.png
diff --git a/markdown/DAE/figures/water_tube_24_1.png b/markdown/DAE/figures/water_tube_24_1.png
diff --git a/markdown/DAE/figures/wheelset_15_1.png b/markdown/DAE/figures/wheelset_15_1.png
diff --git a/markdown/DAE/figures/wheelset_19_1.png b/markdown/DAE/figures/wheelset_19_1.png
diff --git a/markdown/DAE/figures/wheelset_20_1.png b/markdown/DAE/figures/wheelset_20_1.png
diff --git a/markdown/DAE/figures/wheelset_21_1.png b/markdown/DAE/figures/wheelset_21_1.png
diff --git a/markdown/DAE/figures/wheelset_22_1.png b/markdown/DAE/figures/wheelset_22_1.png
diff --git a/markdown/DAE/figures/wheelset_23_1.png b/markdown/DAE/figures/wheelset_23_1.png
diff --git a/markdown/DAE/slider_crank.md b/markdown/DAE/slider_crank.md
diff --git a/markdown/DAE/two_bit_adder.md b/markdown/DAE/two_bit_adder.md
diff --git a/markdown/DAE/water_tube.md b/markdown/DAE/water_tube.md
diff --git a/markdown/DAE/wheelset.md b/markdown/DAE/wheelset.md
diff --git a/script/DAE/ChemicalAkzoNobel.jl b/script/DAE/ChemicalAkzoNobel.jl
@@ -1,19 +1,19 @@
-
 using OrdinaryDiffEq, DiffEqDevTools, Sundials, ModelingToolkit, ODEInterfaceDiffEq,
-      Plots, DASSL, DASKR
+    Plots, DASSL, DASKR
+using OrdinaryDiffEqBDF, OrdinaryDiffEqFIRK, OrdinaryDiffEqRosenbrock, OrdinaryDiffEqSDIRK
 using LinearAlgebra
 using ModelingToolkit: t_nounits as t, D_nounits as D
 
 ModelingToolkit.@parameters begin
-    k₁=18.7
-    k₂=0.58
-    k₃=0.09
-    k₄=0.42
-    kbig=34.4
-    kla=3.3
-    ks=115.83
-    po2=0.9
-    hen=737
+    k₁ = 18.7
+    k₂ = 0.58
+    k₃ = 0.09
+    k₄ = 0.42
+    kbig = 34.4
+    kla = 3.3
+    ks = 115.83
+    po2 = 0.9
+    hen = 737
 end
 
 @variables begin
@@ -22,190 +22,256 @@ end
     y₃(t) = 0.0
     y₄(t) = 0.007
     y₅(t) = 1.0
-    y₆(t) = 115.83*0.444*0.007 # ks*y₁*y₄
+    y₆(t) = 115.83 * 0.444 * 0.007 # ks*y₁*y₄
 end
 
 r₁ = k₁ * (y₁^4.0) * sqrt(abs(y₂))
 r₂ = k₂ * y₃ * y₄
-r₃ = k₂/kbig * y₁ * y₅
-r₄ = k₃*y₁*(y₄^2)
-r₅ = k₄*(y₆^2)*sqrt(abs(y₂))
-fin = kla*(po2/hen-y₂)
-
-eqs = [D(y₁) ~ -2.0 * r₁ + r₂ - r₃ - r₄
-       D(y₂) ~ -0.5 * r₁ - r₄ - 0.5*r₅ + fin
-       D(y₃) ~ r₁ - r₂ + r₃
-       D(y₄) ~ -r₂ + r₃ - 2.0 * r₄
-       D(y₅) ~ r₂ - r₃ + r₅
-       0.0 ~ ks * y₁ * y₄ - y₆]
+r₃ = k₂ / kbig * y₁ * y₅
+r₄ = k₃ * y₁ * (y₄^2)
+r₅ = k₄ * (y₆^2) * sqrt(abs(y₂))
+fin = kla * (po2 / hen - y₂)
+
+eqs = [
+    D(y₁) ~ -2.0 * r₁ + r₂ - r₃ - r₄
+    D(y₂) ~ -0.5 * r₁ - r₄ - 0.5 * r₅ + fin
+    D(y₃) ~ r₁ - r₂ + r₃
+    D(y₄) ~ -r₂ + r₃ - 2.0 * r₄
+    D(y₅) ~ r₂ - r₃ + r₅
+    0.0 ~ ks * y₁ * y₄ - y₆
+]
 
-ModelingToolkit.@mtkbuild sys = ModelingToolkit.ODESystem(eqs, t)
+# @mtkcompile drops y₆ (algebraic). Residual DAE form is hand-written so IDA
+# runs on the original 6-variable index-1 system (MTK DAEProblem + IDA hits a
+# KINSOL handle double-free on finalize for the reduced pure-ODE system).
+ModelingToolkit.@mtkcompile sys = ModelingToolkit.System(eqs, t)
 
 tspan = (0.0, 180.0)
-mtkprob = ODEProblem(sys, [], tspan)
-sol = solve(mtkprob, Rodas4(), abstol = 1/10^14, reltol = 1/10^14)
+mtkprob = ODEProblem(sys, [], tspan; warn_initialize_determined = false)
+sol = solve(mtkprob, Rodas4(), abstol = 1 / 10^14, reltol = 1 / 10^14)
 
-odaeprob = ODAEProblem(sys, [], tspan)
-ode_ref_sol = solve(odaeprob, CVODE_BDF(), abstol = 1/10^14, reltol = 1/10^14);
-
-du = mtkprob.f(mtkprob.u0, mtkprob.p, 0.0)
-du0 = D.(unknowns(sys)) .=> du
-daeprob = DAEProblem(sys, du0, [], tspan)
-ref_sol = solve(daeprob, IDA(), abstol = 1/10^14, reltol = 1/10^14);
+odaeprob = ODEProblem(sys, [], tspan; warn_initialize_determined = false)
+ode_ref_sol = solve(odaeprob, CVODE_BDF(), abstol = 1 / 10^14, reltol = 1 / 10^14)
 
 function akzo(du, u, p, t)
     y₁, y₂, y₃, y₄, y₅, y₆ = u
-    k₁=18.7
-    k₂=0.58
-    k₃=0.09
-    k₄=0.42
-    kbig=34.4
-    kla=3.3
-    ks=115.83
-    po2=0.9
-    hen=737
+    k₁ = 18.7
+    k₂ = 0.58
+    k₃ = 0.09
+    k₄ = 0.42
+    kbig = 34.4
+    kla = 3.3
+    ks = 115.83
+    po2 = 0.9
+    hen = 737
 
     r₁ = k₁ * (y₁^4.0) * sqrt(abs(y₂))
     r₂ = k₂ * y₃ * y₄
-    r₃ = k₂/kbig * y₁ * y₅
-    r₄ = k₃*y₁*(y₄^2)
-    r₅ = k₄*(y₆^2)*sqrt(abs(y₂))
-    fin = kla*(po2/hen-y₂)
+    r₃ = k₂ / kbig * y₁ * y₅
+    r₄ = k₃ * y₁ * (y₄^2)
+    r₅ = k₄ * (y₆^2) * sqrt(abs(y₂))
+    fin = kla * (po2 / hen - y₂)
 
     du[1] = -2.0 * r₁ + r₂ - r₃ - r₄
-    du[2] = -0.5 * r₁ - r₄ - 0.5*r₅ + fin
+    du[2] = -0.5 * r₁ - r₄ - 0.5 * r₅ + fin
     du[3] = r₁ - r₂ + r₃
     du[4] = -r₂ + r₃ - 2.0 * r₄
     du[5] = r₂ - r₃ + r₅
     du[6] = ks * y₁ * y₄ - y₆
-    nothing
+    return nothing
+end
+
+function akzo_dae!(res, du, u, p, t)
+    y₁, y₂, y₃, y₄, y₅, y₆ = u
+    k₁ = 18.7
+    k₂ = 0.58
+    k₃ = 0.09
+    k₄ = 0.42
+    kbig = 34.4
+    kla = 3.3
+    ks = 115.83
+    po2 = 0.9
+    hen = 737
+
+    r₁ = k₁ * (y₁^4.0) * sqrt(abs(y₂))
+    r₂ = k₂ * y₃ * y₄
+    r₃ = k₂ / kbig * y₁ * y₅
+    r₄ = k₃ * y₁ * (y₄^2)
+    r₅ = k₄ * (y₆^2) * sqrt(abs(y₂))
+    fin = kla * (po2 / hen - y₂)
+
+    res[1] = du[1] - (-2.0 * r₁ + r₂ - r₃ - r₄)
+    res[2] = du[2] - (-0.5 * r₁ - r₄ - 0.5 * r₅ + fin)
+    res[3] = du[3] - (r₁ - r₂ + r₃)
+    res[4] = du[4] - (-r₂ + r₃ - 2.0 * r₄)
+    res[5] = du[5] - (r₂ - r₃ + r₅)
+    res[6] = ks * y₁ * y₄ - y₆
+    return nothing
 end
-M = Matrix{Float64}(I, 6, 6);
-M[6, 6] = 0;
+
+u0_akzo = [0.444, 0.00123, 0.0, 0.007, 1.0, 115.83 * 0.444 * 0.007]
+du0_akzo = zeros(6)
+akzo(du0_akzo, u0_akzo, nothing, 0.0)
+du0_akzo[6] = 0.0
+daeprob = DAEProblem(
+    akzo_dae!, du0_akzo, u0_akzo, tspan;
+    differential_vars = [true, true, true, true, true, false]
+)
+ref_sol = solve(daeprob, IDA(), abstol = 1 / 10^14, reltol = 1 / 10^14)
+
+M = Matrix{Float64}(I, 6, 6)
+M[6, 6] = 0
 mmf = ODEFunction(akzo, mass_matrix = M)
-mmprob = ODEProblem(mmf, [0.444, 0.00123, 0.0, 0.007, 1.0, 115.83*0.444*0.007], tspan)
-mm_refsol = solve(mmprob, Rodas5(), reltol = 1e-12, abstol = 1e-12)
+mmprob = ODEProblem(mmf, u0_akzo, tspan)
+mm_refsol = solve(mmprob, Rodas5(), reltol = 1.0e-12, abstol = 1.0e-12)
 
+# mtkprob/odaeprob are the 5-state reduced ODE; daeprob/mmprob are the
+# original 6-state index-1 form. Pair each with a matching reference.
 probs = [mtkprob, daeprob, odaeprob, mmprob]
-refs = [ref_sol, ref_sol, ode_ref_sol, mm_refsol];
+refs = [ode_ref_sol, ref_sol, ode_ref_sol, mm_refsol]
 
 
-plot(ref_sol, idxs = [y₁, y₂, y₃, y₄, y₅, y₆])
+plot(ref_sol)
 
 
 plot(mm_refsol)
 
 
 abstols = 1.0 ./ 10.0 .^ (5:8)
 reltols = 1.0 ./ 10.0 .^ (1:4);
-setups = [Dict(:prob_choice => 1, :alg=>Rosenbrock23()),
-    Dict(:prob_choice => 1, :alg=>Rodas4()),
-    Dict(:prob_choice => 1, :alg=>FBDF()),
-    Dict(:prob_choice => 1, :alg=>QNDF()),
-    Dict(:prob_choice => 1, :alg=>rodas()),
-    Dict(:prob_choice => 1, :alg=>radau()),
-    Dict(:prob_choice => 1, :alg=>RadauIIA5()),
-    Dict(:prob_choice => 2, :alg=>DFBDF()),
-    Dict(:prob_choice => 2, :alg=>IDA())
+setups = [
+    Dict(:prob_choice => 1, :alg => Rosenbrock23()),
+    Dict(:prob_choice => 1, :alg => Rodas4()),
+    Dict(:prob_choice => 1, :alg => FBDF()),
+    Dict(:prob_choice => 1, :alg => QNDF()),
+    Dict(:prob_choice => 1, :alg => NordsieckBDF()),
+    Dict(:prob_choice => 1, :alg => rodas()),
+    Dict(:prob_choice => 1, :alg => radau()),
+    Dict(:prob_choice => 1, :alg => RadauIIA5()),
+    Dict(:prob_choice => 2, :alg => DFBDF()),
+    Dict(:prob_choice => 2, :alg => DNordsieckBDF()),
+    Dict(:prob_choice => 2, :alg => IDA()),
 ]
 
-wp = WorkPrecisionSet(probs, abstols, reltols, setups;
-    save_everystep = false, appxsol = refs, maxiters = Int(1e5), numruns = 10)
+wp = WorkPrecisionSet(
+    probs, abstols, reltols, setups;
+    save_everystep = false, appxsol = refs, maxiters = Int(1.0e5), numruns = 10
+)
 plot(wp)
 
 
 abstols = 1.0 ./ 10.0 .^ (6:8)
 reltols = 1.0 ./ 10.0 .^ (2:4);
-setups = [Dict(:prob_choice => 1, :alg=>Rosenbrock23()),
-    Dict(:prob_choice => 1, :alg=>Rodas4()),
-    Dict(:prob_choice => 2, :alg=>IDA()),
-    Dict(:prob_choice => 3, :alg=>Rosenbrock23()),
-    Dict(:prob_choice => 3, :alg=>Rodas4()),
-    Dict(:prob_choice => 3, :alg=>CVODE_BDF()),
-    Dict(:prob_choice => 3, :alg=>TRBDF2()),
-    Dict(:prob_choice => 3, :alg=>KenCarp4()),
-    Dict(:prob_choice => 4, :alg=>Rodas4())
+setups = [
+    Dict(:prob_choice => 1, :alg => Rosenbrock23()),
+    Dict(:prob_choice => 1, :alg => Rodas4()),
+    Dict(:prob_choice => 2, :alg => IDA()),
+    Dict(:prob_choice => 3, :alg => Rosenbrock23()),
+    Dict(:prob_choice => 3, :alg => Rodas4()),
+    Dict(:prob_choice => 3, :alg => CVODE_BDF()),
+    Dict(:prob_choice => 3, :alg => TRBDF2()),
+    Dict(:prob_choice => 3, :alg => KenCarp4()),
+    Dict(:prob_choice => 4, :alg => Rodas4()),
 ]
-wp = WorkPrecisionSet(probs, abstols, reltols, setups;
-    save_everystep = false, appxsol = refs, maxiters = Int(1e5), numruns = 10)
+wp = WorkPrecisionSet(
+    probs, abstols, reltols, setups;
+    save_everystep = false, appxsol = refs, maxiters = Int(1.0e5), numruns = 10
+)
 plot(wp)
 
 
 abstols = 1.0 ./ 10.0 .^ (6:8)
 reltols = 1.0 ./ 10.0 .^ (3:5);
-setups = [Dict(:prob_choice => 3, :alg=>Rosenbrock23()),
-    Dict(:prob_choice => 3, :alg=>Rodas4()),
-    Dict(:prob_choice => 2, :alg=>IDA()),
-    Dict(:prob_choice => 2, :alg=>DASSL.dassl()),
-    Dict(:prob_choice => 2, :alg=>DASKR.daskr())
+setups = [
+    Dict(:prob_choice => 3, :alg => Rosenbrock23()),
+    Dict(:prob_choice => 3, :alg => Rodas4()),
+    Dict(:prob_choice => 2, :alg => IDA()),
+    Dict(:prob_choice => 2, :alg => DASSL.dassl()),
+    Dict(:prob_choice => 2, :alg => DASKR.daskr()),
 ]
-wp = WorkPrecisionSet(probs, abstols, reltols, setups;
-    save_everystep = false, appxsol = refs, maxiters = Int(1e5), numruns = 10)
+wp = WorkPrecisionSet(
+    probs, abstols, reltols, setups;
+    save_everystep = false, appxsol = refs, maxiters = Int(1.0e5), numruns = 10
+)
 plot(wp)
 
 
 abstols = 1.0 ./ 10.0 .^ (5:8)
 reltols = 1.0 ./ 10.0 .^ (1:4);
-setups = [Dict(:prob_choice => 1, :alg=>Rosenbrock23()),
-    Dict(:prob_choice => 1, :alg=>Rodas4()),
-    Dict(:prob_choice => 1, :alg=>FBDF()),
-    Dict(:prob_choice => 1, :alg=>QNDF()),
-    Dict(:prob_choice => 1, :alg=>rodas()),
-    Dict(:prob_choice => 1, :alg=>radau()),
-    Dict(:prob_choice => 1, :alg=>RadauIIA5()),
-    Dict(:prob_choice => 2, :alg=>DFBDF()),
-    Dict(:prob_choice => 2, :alg=>IDA())
+setups = [
+    Dict(:prob_choice => 1, :alg => Rosenbrock23()),
+    Dict(:prob_choice => 1, :alg => Rodas4()),
+    Dict(:prob_choice => 1, :alg => FBDF()),
+    Dict(:prob_choice => 1, :alg => QNDF()),
+    Dict(:prob_choice => 1, :alg => NordsieckBDF()),
+    Dict(:prob_choice => 1, :alg => rodas()),
+    Dict(:prob_choice => 1, :alg => radau()),
+    Dict(:prob_choice => 1, :alg => RadauIIA5()),
+    Dict(:prob_choice => 2, :alg => DFBDF()),
+    Dict(:prob_choice => 2, :alg => DNordsieckBDF()),
+    Dict(:prob_choice => 2, :alg => IDA()),
 ]
-wp = WorkPrecisionSet(probs, abstols, reltols, setups; error_estimate = :l2,
-    save_everystep = false, appxsol = refs, maxiters = Int(1e5), numruns = 10)
+wp = WorkPrecisionSet(
+    probs, abstols, reltols, setups; error_estimate = :l2,
+    save_everystep = false, appxsol = refs, maxiters = Int(1.0e5), numruns = 10
+)
 plot(wp)
 
 
 abstols = 1.0 ./ 10.0 .^ (6:8)
 reltols = 1.0 ./ 10.0 .^ (2:4);
-setups = [Dict(:prob_choice => 1, :alg=>Rosenbrock23()),
-    Dict(:prob_choice => 1, :alg=>Rodas4()),
-    Dict(:prob_choice => 2, :alg=>IDA()),
-    Dict(:prob_choice => 3, :alg=>Rosenbrock23()),
-    Dict(:prob_choice => 3, :alg=>Rodas4()),
-    Dict(:prob_choice => 3, :alg=>CVODE_BDF()),
-    Dict(:prob_choice => 3, :alg=>TRBDF2()),
-    Dict(:prob_choice => 3, :alg=>KenCarp4())
+setups = [
+    Dict(:prob_choice => 1, :alg => Rosenbrock23()),
+    Dict(:prob_choice => 1, :alg => Rodas4()),
+    Dict(:prob_choice => 2, :alg => IDA()),
+    Dict(:prob_choice => 3, :alg => Rosenbrock23()),
+    Dict(:prob_choice => 3, :alg => Rodas4()),
+    Dict(:prob_choice => 3, :alg => CVODE_BDF()),
+    Dict(:prob_choice => 3, :alg => TRBDF2()),
+    Dict(:prob_choice => 3, :alg => KenCarp4()),
 ]
-wp = WorkPrecisionSet(probs, abstols, reltols, setups; error_estimate = :l2,
-    save_everystep = false, appxsol = refs, maxiters = Int(1e5), numruns = 10)
+wp = WorkPrecisionSet(
+    probs, abstols, reltols, setups; error_estimate = :l2,
+    save_everystep = false, appxsol = refs, maxiters = Int(1.0e5), numruns = 10
+)
 plot(wp)
 
 
 abstols = 1.0 ./ 10.0 .^ (7:12)
 reltols = 1.0 ./ 10.0 .^ (4:9)
 
-setups = [Dict(:prob_choice => 1, :alg=>Rodas5()),
-    Dict(:prob_choice => 3, :alg=>Rodas5()),
-    Dict(:prob_choice => 4, :alg=>Rodas5()),
-    Dict(:prob_choice => 1, :alg=>Rodas4()),
-    Dict(:prob_choice => 3, :alg=>Rodas4()),
-    Dict(:prob_choice => 4, :alg=>Rodas4()),
-    Dict(:prob_choice => 1, :alg=>FBDF()),
-    Dict(:prob_choice => 1, :alg=>QNDF()),
-    Dict(:prob_choice => 1, :alg=>rodas()),
-    Dict(:prob_choice => 1, :alg=>radau()),
-    Dict(:prob_choice => 1, :alg=>RadauIIA5()),
-    Dict(:prob_choice => 2, :alg=>DFBDF()),
-    Dict(:prob_choice => 2, :alg=>IDA()),
-    Dict(:prob_choice => 2, :alg=>DASKR.daskr())
+setups = [
+    Dict(:prob_choice => 1, :alg => Rodas5()),
+    Dict(:prob_choice => 3, :alg => Rodas5()),
+    Dict(:prob_choice => 4, :alg => Rodas5()),
+    Dict(:prob_choice => 1, :alg => Rodas4()),
+    Dict(:prob_choice => 3, :alg => Rodas4()),
+    Dict(:prob_choice => 4, :alg => Rodas4()),
+    Dict(:prob_choice => 1, :alg => FBDF()),
+    Dict(:prob_choice => 1, :alg => QNDF()),
+    Dict(:prob_choice => 1, :alg => NordsieckBDF()),
+    Dict(:prob_choice => 1, :alg => rodas()),
+    Dict(:prob_choice => 1, :alg => radau()),
+    Dict(:prob_choice => 1, :alg => RadauIIA5()),
+    Dict(:prob_choice => 2, :alg => DFBDF()),
+    Dict(:prob_choice => 2, :alg => DNordsieckBDF()),
+    Dict(:prob_choice => 2, :alg => IDA()),
+    Dict(:prob_choice => 2, :alg => DASKR.daskr()),
 ]
 
-wp = WorkPrecisionSet(probs, abstols, reltols, setups;
-    save_everystep = false, appxsol = refs, maxiters = Int(1e5), numruns = 10)
+wp = WorkPrecisionSet(
+    probs, abstols, reltols, setups;
+    save_everystep = false, appxsol = refs, maxiters = Int(1.0e5), numruns = 10
+)
 plot(wp)
 
 
-wp = WorkPrecisionSet(probs, abstols, reltols, setups; error_estimate = :l2,
-    save_everystep = false, appxsol = refs, maxiters = Int(1e5), numruns = 10)
+wp = WorkPrecisionSet(
+    probs, abstols, reltols, setups; error_estimate = :l2,
+    save_everystep = false, appxsol = refs, maxiters = Int(1.0e5), numruns = 10
+)
 plot(wp)
 
 
 using SciMLBenchmarks
 SciMLBenchmarks.bench_footer(WEAVE_ARGS[:folder], WEAVE_ARGS[:file])
-