function dist(ds::CDS, d0, Ttr, T; diff_eq_kwargs=Dict(:abstol=>1e-14,

    :reltol=>1e-14, :solver=>Vern9(), :maxiters=>1e9),

    inittest=ChaosTools.inittest_default(dimension(ds)))



    pinteg = parallel_integrator(ds,

            [deepcopy(get_state(ds)), inittest(get_state(ds), d0)];

            diff_eq_kwargs=diff_eq_kwargs)



    if Ttr != 0

        add_tstop!(pinteg, Ttr)

        pinteg.opts.stop_at_next_tstop = true

        for (u,t) in tuples(pinteg) end

    end



    save_func(u, t, integrator) = norm(u[1] - u[2])

    saved_values = SavedValues(eltype(ds.prob.tspan), eltype(ds.prob.u0[1]))

    cb = SavingCallback(save_func, saved_values)



    prob = remake(pinteg.sol.prob, u0=pinteg.u, tspan=(Ttr, Ttr + T), callback=cb)

    solver, newkw = DynamicalSystemsBase.extract_solver(diff_eq_kwargs)

    sol = solve(prob, solver; newkw...)



    return saved_values.saveval, saved_values.t

end