Restore Shampine Rosenbrock23/32 stiff interpolation

Replace H=zeros(T,0,3) (Hermite fallback) with the 1x3 dense coefficient
matrix H=[0, -1/(gamma*(1-2*gamma)), 0] for both Rosenbrock23 and Rosenbrock32,
yielding interp_order=1 and one dense stage vector.

Add interp_order==1 branches to all eight _ode_interpolant/_ode_interpolant!
overloads (Val{0}/Val{1}, idxs::Nothing/idxs, returning/mutating), which
implement p(Theta)=(1-Theta)*y0+Theta*y1+Theta*(1-Theta)*K[1] and its
first derivative.

This algebraically recovers the original MATLAB ODE Suite Shampine
'free' 2nd-order stiff-aware dense output. Fixes regression noticed
by ChrisRackauckas in SciML#3273.

Author: Harsh Singh

lib/OrdinaryDiffEqRosenbrock/src/rosenbrock_interpolants.jl

@@ -40,6 +40,9 @@ end
         # Generic Hermite: k[1]=f₀, k[2]=f₁ (methods with no H matrix)
         @.. Θ1 * y₀ + Θ * y₁ +
             Θ * (Θ - 1) * ((1 - 2Θ) * (y₁ - y₀) + (Θ - 1) * dt * k[1] + Θ * dt * k[2])
+    elseif hasproperty(cache, :interp_order) && cache.interp_order == 1
+        # Shampine Rosenbrock23/32: p(Θ) = (1-Θ)y₀ + Θy₁ + Θ(1-Θ)K[1]
+        @.. Θ1 * y₀ + Θ * (y₁ + Θ1 * k[1])
     elseif !hasproperty(cache, :interp_order) || cache.interp_order == 2
         @.. Θ1 * y₀ + Θ * (y₁ + Θ1 * (k[1] + Θ * k[2]))
     elseif cache.interp_order == 4
@@ -62,6 +65,8 @@ end
     if hasproperty(cache, :interp_order) && cache.interp_order == -1
         @views @.. Θ1 * y₀[idxs] + Θ * y₁[idxs] +
             Θ * (Θ - 1) * ((1 - 2Θ) * (y₁[idxs] - y₀[idxs]) + (Θ - 1) * dt * k[1][idxs] + Θ * dt * k[2][idxs])
+    elseif hasproperty(cache, :interp_order) && cache.interp_order == 1
+        @views @.. Θ1 * y₀[idxs] + Θ * (y₁[idxs] + Θ1 * k[1][idxs])
     elseif !hasproperty(cache, :interp_order) || cache.interp_order == 2
         @views @.. Θ1 * y₀[idxs] + Θ * (y₁[idxs] + Θ1 * (k[1][idxs] + Θ * k[2][idxs]))
     elseif cache.interp_order == 4
@@ -90,6 +95,8 @@ end
     if hasproperty(cache, :interp_order) && cache.interp_order == -1
         @.. out = Θ1 * y₀ + Θ * y₁ +
             Θ * (Θ - 1) * ((1 - 2Θ) * (y₁ - y₀) + (Θ - 1) * dt * k[1] + Θ * dt * k[2])
+    elseif hasproperty(cache, :interp_order) && cache.interp_order == 1
+        @.. out = Θ1 * y₀ + Θ * (y₁ + Θ1 * k[1])
     elseif !hasproperty(cache, :interp_order) || cache.interp_order == 2
         @.. out = Θ1 * y₀ + Θ * (y₁ + Θ1 * (k[1] + Θ * k[2]))
     elseif cache.interp_order == 4
@@ -113,6 +120,8 @@ end
     if hasproperty(cache, :interp_order) && cache.interp_order == -1
         @views @.. out = Θ1 * y₀[idxs] + Θ * y₁[idxs] +
             Θ * (Θ - 1) * ((1 - 2Θ) * (y₁[idxs] - y₀[idxs]) + (Θ - 1) * dt * k[1][idxs] + Θ * dt * k[2][idxs])
+    elseif hasproperty(cache, :interp_order) && cache.interp_order == 1
+        @views @.. out = Θ1 * y₀[idxs] + Θ * (y₁[idxs] + Θ1 * k[1][idxs])
     elseif !hasproperty(cache, :interp_order) || cache.interp_order == 2
         @views @.. out = Θ1 * y₀[idxs] + Θ * (y₁[idxs] + Θ1 * (k[1][idxs] + Θ * k[2][idxs]))
     elseif cache.interp_order == 4
@@ -151,6 +160,9 @@ end
         @.. 6 * Θ * (1 - Θ) * (y₁ - y₀) / dt +
             k[1] * (3 * Θ^2 - 4 * Θ + 1) +
             k[2] * (3 * Θ^2 - 2 * Θ)
+    elseif hasproperty(cache, :interp_order) && cache.interp_order == 1
+        # Shampine Rosenbrock23/32: dp/dt = (y₁-y₀ + (1-2Θ)K[1]) / dt
+        @.. ((1 - 2 * Θ) * k[1] - y₀ + y₁) / dt
     elseif !hasproperty(cache, :interp_order) || cache.interp_order == 2
         @.. (k[1] + Θ * (-2 * k[1] + 2 * k[2] - 3 * k[2] * Θ) - y₀ + y₁) / dt
     elseif cache.interp_order == 4
@@ -188,6 +200,8 @@ end
         @views @.. 6 * Θ * (1 - Θ) * (y₁[idxs] - y₀[idxs]) / dt +
             k[1][idxs] * (3 * Θ^2 - 4 * Θ + 1) +
             k[2][idxs] * (3 * Θ^2 - 2 * Θ)
+    elseif hasproperty(cache, :interp_order) && cache.interp_order == 1
+        @views @.. ((1 - 2 * Θ) * k[1][idxs] - y₀[idxs] + y₁[idxs]) / dt
     elseif !hasproperty(cache, :interp_order) || cache.interp_order == 2
         @views @.. (
             k[1][idxs] +
@@ -231,6 +245,8 @@ end
         @.. out = 6 * Θ * (1 - Θ) * (y₁ - y₀) / dt +
             k[1] * (3 * Θ^2 - 4 * Θ + 1) +
             k[2] * (3 * Θ^2 - 2 * Θ)
+    elseif hasproperty(cache, :interp_order) && cache.interp_order == 1
+        @.. out = ((1 - 2 * Θ) * k[1] - y₀ + y₁) / dt
     elseif !hasproperty(cache, :interp_order) || cache.interp_order == 2
         @.. out = (k[1] + Θ * (-2 * k[1] + 2 * k[2] - 3 * k[2] * Θ) - y₀ + y₁) / dt
     elseif cache.interp_order == 4
@@ -270,6 +286,8 @@ end
         @views @.. out = 6 * Θ * (1 - Θ) * (y₁[idxs] - y₀[idxs]) / dt +
             k[1][idxs] * (3 * Θ^2 - 4 * Θ + 1) +
             k[2][idxs] * (3 * Θ^2 - 2 * Θ)
+    elseif hasproperty(cache, :interp_order) && cache.interp_order == 1
+        @views @.. out = ((1 - 2 * Θ) * k[1][idxs] - y₀[idxs] + y₁[idxs]) / dt
     elseif !hasproperty(cache, :interp_order) || cache.interp_order == 2
         @views @.. out = (
             k[1][idxs] +

lib/OrdinaryDiffEqRosenbrock/src/rosenbrock_tableaus.jl

@@ -38,7 +38,11 @@ function Rosenbrock23RodasTableau(T, T2)
         convert(T, igamma / 6),
     ]
 
-    H = zeros(T, 0, 3)
+    # Shampine 2nd-order "free" stiff interpolation (Rodas-encoded).
+    # K[1] = H[1,2]*ks[2] = -ks[2]/(γ(1-2γ)); combined with interp_order=1 formula
+    # p(Θ) = (1-Θ)y₀ + Θy₁ + Θ(1-Θ)K[1] this exactly recovers the MATLAB ODE Suite
+    # dense output c1*k̃₁ + c2*k̃₂ with c1=Θ(1-Θ)/(1-2γ), c2=Θ(Θ-2γ)/(1-2γ).
+    H = T[zero(T) convert(T, -1 / (gamma * (1 - 2 * gamma))) zero(T)]
 
     return RodasTableau(A, C, gamma, c, d, H, b, btilde)
 end
@@ -75,7 +79,8 @@ function Rosenbrock32RodasTableau(T, T2)
         convert(T, igamma / 6),
     ]
 
-    H = zeros(T, 0, 3)
+    # Same Shampine interpolation as Rosenbrock23 (identical γ, same stencil)
+    H = T[zero(T) convert(T, -1 / (gamma * (1 - 2 * gamma))) zero(T)]
 
     return RodasTableau(A, C, gamma, c, d, H, b, btilde)
 end