stability_contour.jl

using DiffEqDevTools, BifurcationKit, Setfield, Plots

mutable struct TermCache
    start
    sign::Bool
    count::Int
end

function stability_curve(tab)
    function curve((x, ), (y, ))
        sol = stability_region(x+y*im, tab)
        [abs(sol) - 1]
    end
    cache = TermCache(nothing, true, 0)
    function terminate(z, tau, step, contResult)
        if step == 1
            cache.start = z.p
        elseif step == 2
            cache.sign = cache.start - z.p > 0
        else
            if (cache.start - z.p > 0) != cache.sign
                cache.sign = !cache.sign
                cache.count += 1
            end
            return cache.count <= 1
        end
        return true
    end
    sol, _ = BifurcationKit.continuation(curve, [-0.1,], [0.0], (@lens _[1]), ContinuationPar(pMin=-80.0, pMax=80.0, maxSteps=1000); finaliseSolution = terminate)
    return sol
end
sol_tsit5 = stability_curve(constructTsitouras5())
plot(sol_tsit5.branch[2, :], sol_tsit5.branch[1, :])