def METHOD_NAME(integrate_func, dt):
    p = HernquistPotential(m=1e11, c=0.5, units=galactic)
    H = Hamiltonian(potential=p)
    w0 = np.array(
        [
            [0.0, 10.0, 0.0, 0.2, 0.0, 0.0],
            [10.0, 0.0, 0.0, 0.0, 0.2, 0.0],
            [0.0, 10.0, 0.0, 0.0, 0.0, 0.2],
        ]
    )
    # 1024 steps
    t = np.linspace(0, dt * 1024, 1024 + 1)
    t_all, w_all = integrate_func(H, w0, t)
    t_f, w_f = integrate_func(H, w0, t, store_all=False)
    assert t_all[-1] == t_f[0]
    assert np.allclose(w_all[-1], w_f)