-- dop853.lua: Dormand-Prince 8th order (13 NFE, fixed step)

solver = {
    name        = "dop853",
    display     = "DOP853 (13 NFE)",
    description = "8th order Dormand-Prince (maximum accuracy)",
    nfe         = 13,
    order       = 8,
    needs_model = true,
    stateful    = false,
    stochastic  = false,
}

local C = {
    1/18, 1/12, 1/8, 5/16, 3/8,
    59/400, 93/200, 5490023248/9719169821,
    13/20, 1201146811/1299019798,
    1, 1,
}
local A = {
    {1/18},
    {1/48, 1/16},
    {1/32, 0, 3/32},
    {5/16, 0, -75/64, 75/64},
    {3/80, 0, 0, 3/16, 3/20},
    {29443841/614563906, 0, 0, 77736538/692538347, -28693883/1125000000, 23124283/1800000000},
    {16016141/946692911, 0, 0, 61564180/158732637, 22789713/633445777, 545815736/2771057229, -180193667/1043307555},
    {39632708/573591083, 0, 0, -433636366/683701615, -421739975/2616292301, 100302831/723423059, 790204164/839813087, 800635310/3783071287},
    {246121993/1340847787, 0, 0, -37695042795/15268766246, -309121744/1061227803, -12992083/490766935, 6005943493/2108947869, 393006217/1396673457, 123872331/1001029789},
    {-1028468189/846180014, 0, 0, 8478235783/508512852, 1311729495/1432422823, -10304129995/1701304382, -48777925059/3047939560, 15336726248/1032824649, -45442868181/3398467696, 3065993473/597172653},
    {185892177/718116043, 0, 0, -3185094517/667107341, -477755414/1098053517, -703635378/230739211, 5731566787/1027545527, 5232866602/850066563, -4093664535/808688257, 3962137247/1805957418, 65686358/487910083},
    {403863854/491063109, 0, 0, -5068492393/434740067, -411421997/543043805, 652783627/914296604, 11173962825/925320556, -13158990841/6184727034, 3936647629/1978049680, -160528059/685178525, 248638103/1413531060, 0},
}
local B = {
    14005451/335480064, 0, 0, 0, 0,
    -59238493/1068277825, 181606767/758867731, 561292985/797845732,
    -1041891430/1371343529, 760417239/1151165299, 118820643/751138087,
    -528747749/2220607170, 1/4,
}

function step(xt, vt, t_curr, t_prev, n, model_fn, vt_buf)
    local dt = t_curr - t_prev
    local xt_orig = {}
    local k1 = {}
    for i = 0, n-1 do xt_orig[i] = xt[i]; k1[i] = vt[i] end

    local ks = {k1}
    local x_tmp = {}

    -- 12 extra stages
    for s = 1, 12 do
        local a_row = A[s]
        for i = 0, n-1 do
            local combo = 0
            for j = 1, #a_row do
                if a_row[j] ~= 0 and ks[j] then
                    combo = combo + a_row[j] * ks[j][i]
                end
            end
            x_tmp[i] = xt_orig[i] - dt * combo
        end
        -- Write x_tmp to xt for model_fn
        for i = 0, n-1 do xt[i] = x_tmp[i] end
        model_fn(xt, t_curr - C[s] * dt)
        ks[s+1] = {}
        for i = 0, n-1 do ks[s+1][i] = vt_buf[i] end
    end

    -- Combine
    for i = 0, n-1 do
        local sol = 0
        for j = 1, 13 do
            if B[j] ~= 0 and ks[j] then sol = sol + B[j] * ks[j][i] end
        end
        xt[i] = xt_orig[i] - dt * sol
    end
end