# Verner 9 Runge Kutta Interpolator Issue

**URL:** <https://discourse.julialang.org/t/verner-9-runge-kutta-interpolator-issue/131591>\
**Category:** New to Julia\
**Created:** [August 13, 2025, 5:26pm UTC](https://discourse.julialang.org/t/verner-9-runge-kutta-interpolator-issue/131591 "2025-08-13T17:26:46Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![may7even](https://avatars.discourse-cdn.com/v4/letter/m/8dc957/32.png) [@may7even](https://discourse.julialang.org/u/may7even)\
**Post date:** [August 13, 2025, 5:26pm UTC](https://discourse.julialang.org/t/verner-9-runge-kutta-interpolator-issue/131591/1 "2025-08-13T17:26:46Z")

</div>

Hello, I am trying to implement the vern9 runge kutta method for a side project, and so far, everything has been good, except for the interpolation of the states for a dense output. For some reason, when compared with a true data, the interpolated values are off by quite a lot. I am wondering if my understanding of how to implement the interpolator is just wrong then. It seems straightforward from the source code online, and also from the vern 1993 paper as well. Any help would be greatly appreciated!

```julia-auto
def vern9pre0(theta, r):

    theta2 = theta * theta
    
    b1 = theta * evalpoly(theta, r["r011"], r["r012"], r["r013"], r["r014"], r["r015"], r["r016"], r["r017"], r["r018"], r["r019"])
    b8 = theta2 * evalpoly(theta, r["r082"], r["r083"], r["r084"], r["r085"], r["r086"], r["r087"], r["r088"], r["r089"])
    b9 = theta2 * evalpoly(theta, r["r092"], r["r093"], r["r094"], r["r095"], r["r096"], r["r097"], r["r098"], r["r099"])
    b10 = theta2 * evalpoly(theta, r["r102"], r["r103"], r["r104"], r["r105"], r["r106"], r["r107"], r["r108"], r["r109"])
    b11 = theta2 * evalpoly(theta, r["r112"], r["r113"], r["r114"], r["r115"], r["r116"], r["r117"], r["r118"], r["r119"])
    b12 = theta2 * evalpoly(theta, r["r122"], r["r123"], r["r124"], r["r125"], r["r126"], r["r127"], r["r128"], r["r129"])
    b13 = theta2 * evalpoly(theta, r["r132"], r["r133"], r["r134"], r["r135"], r["r136"], r["r137"], r["r138"], r["r139"])
    b14 = theta2 * evalpoly(theta, r["r142"], r["r143"], r["r144"], r["r145"], r["r146"], r["r147"], r["r148"], r["r149"])
    b15 = theta2 * evalpoly(theta, r["r152"], r["r153"], r["r154"], r["r155"], r["r156"], r["r157"], r["r158"], r["r159"])
    b17 = theta2 * evalpoly(theta, r["r172"], r["r173"], r["r174"], r["r175"], r["r176"], r["r177"], r["r178"], r["r179"])
    b18 = theta2 * evalpoly(theta, r["r182"], r["r183"], r["r184"], r["r185"], r["r186"], r["r187"], r["r188"], r["r189"])
    b19 = theta2 * evalpoly(theta, r["r192"], r["r193"], r["r194"], r["r195"], r["r196"], r["r197"], r["r198"], r["r199"])
    b20 = theta2 * evalpoly(theta, r["r202"], r["r203"], r["r204"], r["r205"], r["r206"], r["r207"], r["r208"], r["r209"])
    b21 = theta2 * evalpoly(theta, r["r212"], r["r213"], r["r214"], r["r215"], r["r216"], r["r217"], r["r218"], r["r219"])
    b22 = theta2 * evalpoly(theta, r["r222"], r["r223"], r["r224"], r["r225"], r["r226"], r["r227"], r["r228"], r["r229"])
    b23 = theta2 * evalpoly(theta, r["r232"], r["r233"], r["r234"], r["r235"], r["r236"], r["r237"], r["r238"], r["r239"])
    b24 = theta2 * evalpoly(theta, r["r242"], r["r243"], r["r244"], r["r245"], r["r246"], r["r247"], r["r248"], r["r249"])
    b25 = theta2 * evalpoly(theta, r["r252"], r["r253"], r["r254"], r["r255"], r["r256"], r["r257"], r["r258"], r["r259"])
    b26 = theta2 * evalpoly(theta, r["r262"], r["r263"], r["r264"], r["r265"], r["r266"], r["r267"], r["r268"], r["r269"])
    
    return (b1, b8, b9, b10, b11, b12, b13, b14, b15,
            b17, b18, b19, b20, b21, b22, b23, b24, b25, b26)

def ode_interpolant(theta, dt, y0, y1, k, r):
    (b1, b8, b9, b10, b11, b12, b13, b14, b15,
     b17, b18, b19, b20, b21, b22, b23, b24, b25, b26) = vern9pre0(theta, r)

    return y0 + dt * (
        k[0] * b1 + k[1] * b8 + k[2] * b9 + k[3] * b10 +
        k[4] * b11 + k[5] * b12 + k[6] * b13 + k[7] * b14 +
        k[8] * b15 + k[10] * b17 + k[11] * b18 + k[12] * b19 +
        k[13] * b20 + k[14] * b21 + k[15] * b22 + k[16] * b23 +
        k[17] * b24 + k[18] * b25 + k[19] * b26
    )

```

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [August 13, 2025, 6:02pm UTC](https://discourse.julialang.org/t/verner-9-runge-kutta-interpolator-issue/131591/2 "2025-08-13T18:02:11Z")

</div>

Are you sure those are the right k’s? It won’t be the same k’s as the integration.

---

<div class="post-metadata">

**Author:** ![may7even](https://avatars.discourse-cdn.com/v4/letter/m/8dc957/32.png) [@may7even](https://discourse.julialang.org/u/may7even)\
**Post date:** [August 13, 2025, 7:22pm UTC](https://discourse.julialang.org/t/verner-9-runge-kutta-interpolator-issue/131591/3 "2025-08-13T19:22:30Z")

</div>

Yeah, the k values are mapped according to what’s in the source code. I’m pretty sure the mistake is on my end, but I wanted to be sure of my understanding of the theory.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [August 14, 2025, 4:40am UTC](https://discourse.julialang.org/t/verner-9-runge-kutta-interpolator-issue/131591/4 "2025-08-14T04:40:09Z")

</div>

So yeah if you mapped the `k`s correctly to their right indexes then this is at least right in theory, sans any typos I didn’t catch.