You can use normalize_position
to wrap any evaluated position within your model boundaries, e.g. in the line you commented out in your agent_step!
function you could write instead
newpos = normalize_position(agent.pos + dt*agent.desired_velocity, model)
move_agent!(agent, newpos, model)
So you can just use whatever integrator you want, just call normalize_position
before trying to update the agent position.
Alternatively, you can use the walk!
function instead which respects the model boundary conditions without the need for you to do it explicitly.