How I was originally doing https://github.com/alecmce/dust/blob/physics/src/dust/physics/systems/PhysicsSystem.hx. It worked well enough, but had gradual errors that I attributed to deriving A(t+dt) based on B(t) rather than based on B(t+dt).

RK4 Integrator First Try

package dust.altphysics.systems;

import dust.altphysics.data.Physics;
import dust.altphysics.data.Physics;
import dust.systems.System;
import dust.altphysics.data.State;
import dust.altphysics.data.Physics;
import dust.altphysics.data.Derivative;
import dust.collections.api.Collection;

class PhysicsSystem implements System
{
    @inject public var collection:Collection;

    inline static var ONE_SIXTH:Float = 1.0 / 6.0;

    public function new() {}

    public function start() {}
    public function stop() {}

    public function iterate(deltaTime:Float)
    {
        for (entity in collection)
        {
            var state:State = entity.get(State);
            var physics:Physics = entity.get(Physics);
            applyRK4IntegrationToStateUnderForces(physics, state, deltaTime);
        }
    }

        inline function applyRK4IntegrationToStateUnderForces(physics:Physics, state:State, deltaTime:Float)
        {
            var halfDeltaTime = deltaTime * 0.5;

            deriveChangeOfStateByForces(physics, state, physics.a);
            applyDerivativeToWorkingState(physics, state, physics.a, halfDeltaTime);
            deriveChangeOfStateByForces(physics, physics.workingState, physics.b);
            applyDerivativeToWorkingState(physics, state, physics.b, halfDeltaTime);
            deriveChangeOfStateByForces(physics, physics.workingState, physics.c);
            applyDerivativeToWorkingState(physics, state, physics.c, deltaTime);
            deriveChangeOfStateByForces(physics, physics.workingState, physics.d);

            updateStateByInterpolatingDerivatives(physics, state, deltaTime);
        }

            inline function deriveChangeOfStateByForces(physics:Physics, state:State, derivative:Derivative)
            {
                derivative.forceX = 0.0;
                derivative.forceY = 0.0;
                derivative.torque = 0.0;

                derivative.velocityX = state.velocityX;
                derivative.velocityY = state.velocityY;
                derivative.spin = state.spin;

                physics.apply(state, derivative);
            }

            inline function applyDerivativeToWorkingState(physics:Physics, state:State, derivative:Derivative, deltaTime:Float)
            {
                var workingState = physics.workingState;
                workingState.read(state);
                workingState.positionX += derivative.velocityX * deltaTime;
                workingState.positionY += derivative.velocityY * deltaTime;
                workingState.rotation += derivative.spin * deltaTime;

                workingState.momentumX += derivative.forceX * deltaTime;
                workingState.momentumY += derivative.forceY * deltaTime;
                workingState.angularMomentum += derivative.torque * deltaTime;

                workingState.updateVelocityFromMomentum();
            }

            inline function updateStateByInterpolatingDerivatives(physics:Physics, state:State, deltaTime:Float)
            {
                var a = physics.a;
                var b = physics.b;
                var c = physics.c;
                var d = physics.d;

                var sixthTime = deltaTime * ONE_SIXTH;

                state.positionX += sixthTime * (a.velocityX + 2.0 * (b.velocityX + c.velocityX) + d.velocityX);
                state.positionY += sixthTime * (a.velocityY + 2.0 * (b.velocityY + c.velocityY) + d.velocityY);
                state.rotation += sixthTime * (a.spin + 2.0 * (b.spin + c.spin) + d.spin);

                state.momentumX += sixthTime * (a.forceX + 2.0 * (b.forceX + c.forceX) + d.forceX);
                state.momentumY += sixthTime * (a.forceY + 2.0 * (b.forceY + c.forceY) + d.forceY);
                state.angularMomentum += sixthTime * (a.torque + 2.0 * (b.torque + c.torque) + d.torque);
                state.updateVelocityFromMomentum();
            }
}