why do you repeat omega? It seems there are mistakes in the parameters. I would suggest using named tuples.
Also, for the period to be 2pi, you need omega = 1.
As a good test for what you are doing, you should plot the residual:
pb(zᵧ, [1.00, 2.5, 0.2, 1.0, 1.00])[1:end-1] |> plot
It shows that the last two equations, for which you have analytical expression, are not behaving as intended. The guess is wrong (4th component) and it should be ordered in the other way (this is not well documented, so … thank you for spotting!). Basically:
zᵧ = vcat(vec(hcat(cos.(γ), -sin.(γ), cos.(γ), sin.(γ))'), 2pi)
To get the solution, it is easier to call:
sol = BifurcationKit.extractTimeSlices(pb, output)
plot(sol[3,:])
You could also use a shooting method, maybe easier to set up from an ODEProblem.