It looks like find_closed_orbit
and newton!
are mishandling the solution after correctly converging:
Specifically, newton!( f!, y, x, ...)
places the final solution into x
, and improperly labels y
as the solution vector. In fact y
is just the residual.
Then, find_closed_orbit(ring; v0=zeros(6))
basically throws away v0 == x
and uses v == y
as the solution.
The fix is to simply replace x <==> y
in the order of arguments of newton!
and then also update its return statements to give (; u=x, ...)