def METHOD_NAME(solvers):
    testsolver, scipysolver, ode = solvers
    # perform steps of the same size
    testsolver.initialize(ode)
    scipysolver.step()
    teststate = diffeq.ODESolverState(
        ivp=ode,
        rv=randvars.Constant(scipysolver.y_old),
        t=scipysolver.t_old,
        error_estimate=None,
        reference_state=None,
    )
    state = testsolver.attempt_step(
        state=teststate, dt=scipysolver.t - scipysolver.t_old
    )
    # sanity check: the steps are the same
    # (this is contained in a different test already, but if this one
    # does not work, the dense output test below is meaningless)
    np.testing.assert_allclose(scipysolver.y, state.rv.mean)
    testsolver_dense = testsolver.dense_output()
    scipy_dense = scipysolver._dense_output_impl()
    t_old = scipysolver.t_old
    t = scipysolver.t
    t_mid = (t_old + t) / 2.0
    for time in [t_old, t, t_mid]:
        test_dense = testsolver_dense(time)
        ref_dense = scipy_dense(time)
        np.testing.assert_allclose(
            test_dense,
            ref_dense,
            atol=1e-13,
            rtol=1e-13,
        )