Hi, This is a great package for bifurcation, thank you.

I am trying to make it work for a simple ODE system, but my hopf continuation failed. Below is my code:

```
function Fitz(u,p)
x,y = u
@unpack I,a,b,c = p
du = similar(u)
du[1] = x-x^3/3-y + I
du[2] = a*(b*x-c*y)
return du
end
sol0 = [0.1,0.1]
par_Fitz = (I=-2.0,a=0.08,b=1.,c=0.8)
ls = GMRESKrylovKit(dim = 100)
optnewton = NewtonPar(tol = 1e-11, verbose = true, linsolver = ls)
out, _, _ = @time newton( Fitz, sol0, par_Fitz, optnewton)
optnew = NewtonPar(verbose = true, tol = 1e-10)
optcont = ContinuationPar(dsmin = 0.01, dsmax = 0.1, ds= 0.05, pMax = 5.1, pMin = -2.1,
newtonOptions = setproperties(optnew; maxIter = 50, tol = 1e-10),
maxSteps = 300, plotEveryStep = 40,
detectBifurcation = 3, nInversion = 4, tolBisectionEigenvalue = 1e-12, dsminBisection = 1e-5)
br, _ = @time continuation(Fitz, out, par_Fitz, (@lens _.I),
optcont; plot = true, verbosity = 0,
printSolution = (x, p) -> x[1])
Ddt(f, x, p, dx) = ForwardDiff.derivative(t -> f(x .+ t .* dx, p), 0.)
dFbpSecBif(x,p) = ForwardDiff.jacobian( z-> Fitz(z,p), x)
d1FbpSecBif(x,p,dx1) = Ddt((z, p0) -> Fitz(z, p0), x, p, dx1)
d2FbpSecBif(x,p,dx1,dx2) = Ddt((z, p0) -> d1FbpSecBif(z, p0, dx1), x, p, dx2)
d3FbpSecBif(x,p,dx1,dx2,dx3) = Ddt((z, p0) -> d2FbpSecBif(z, p0, dx1, dx2), x, p, dx3)
jet = (Fitz, dFbpSecBif, d2FbpSecBif, d3FbpSecBif)
# index of the Hopf point in br.bifpoint
norminf = x -> norm(x, Inf)
ind_hopf = 1
hopfpoint, _, flag = @time newton(Fitz, dFbpSecBif,
br, ind_hopf, par_Fitz, (@lens _.I); normN = norminf)
flag && printstyled(color=:red, "--> We found a Hopf Point at μ = ", hopfpoint.p[1],
", from μ = ", br.bifpoint[ind_hopf].param, "\n")
# automatic branch switching from Hopf point
opt_po = NewtonPar(tol = 1e-5, verbose = true, maxIter = 15)
opts_po_cont = ContinuationPar(dsmin = 0.01, dsmax = 0.1, ds = 0.1, pMax = 5.5,pMin=-2.1, maxSteps = 200,
newtonOptions = opt_po, saveSolEveryStep = 2, plotEveryStep = 1, nev = 11, precisionStability = 1e-4,
detectBifurcation = 3, dsminBisection = 1e-4, maxBisectionSteps = 10);
# number of time slices for the periodic orbit
M = 51
probFD = PeriodicOrbitTrapProblem(M = M)
br_po, = continuation(
# arguments for branch switching from the first
# Hopf bifurcation point
jet..., br, 1,
# arguments for continuation
opts_po_cont, probFD;
# OPTIONAL parameters
# we want to jump on the new branch at phopf + δp
# ampfactor is a factor to increase the amplitude of the guess
δp = 0.05, ampfactor = 1.1,
# specific method for solving linear system
# of Periodic orbits with trapeze method
linearPO = :FullLU,
# regular options for continuation
verbosity = 3, plot = true,
printSolution = (x,p)->x[1], normC = norminf)
```

The branch br is correct, while br_po is not correct. Could you tell what is wrong?