The inconsistent results with the same inputs were due to the adaptive dt in the solver. This issue is solved by resetting the integrator time and dt before stepping:
function f!(xnext, x, u, _, p)
(integrator, set_x, set_u, _) = p
set_t!(integrator, 0.0)
set_proposed_dt!(integrator, Ts)
set_x(integrator, x)
set_u(integrator, u)
step!(integrator, Ts)
xnext .= integrator.u # sol.u is the state, called x in the function
nothing
end