Orbitdiagram DynamicalSystems.jl

Parameter that changes at high values leads to synchronization, that is, a regular mode appears. Ideally, the diagram should show the transition from chaos to regularity. The diagram show this, but there are windows that look as if the trajectory jumps to another attractor (the system has multistability). Is there any way to make sure that inheritance is used ?

Diagram:


Code:

function sigma(x)
    return @fastmath 1.0 / ( 1.0 + exp( -10.0 * ( x  - ( - 0.25 ) ) ) )
end

function HR(u, p, t)
        
    a, b, c, d, s, xr, r,  I, vs, k1, k2, el_link  = p
    x1, y1, z1, x2, y2, z2 = u
    
    du1 = y1 + b * x1 ^ 2 - a * x1 ^3 - z1 + I - k1 * ( x1 - vs ) * sigma(x2) + el_link * ( x2 - x1 )
    du2 = c - d * x1 ^2 - y1
    du3 = r * ( s * ( x1 - xr ) - z1 )
    
    du4 = y2 + b * x2 ^ 2 - a * x2 ^3 - z2 + I - k2 * ( x2 - vs ) * sigma(x1) + el_link * ( x1 - x2 )
    du5 = c - d * x2 ^2 - y2
    du6 = r * ( s * ( x2 - xr ) - z2 )
    return SVector(du1, du2, du3,
                    du4, du5, du6)
end

a = 1.; b = 3.; c = 1.; d = 5.;
xr = -1.6; r = 0.01; s = 5.; I = 4.; xv = 2.;
k1= -0.17; k2 = -0.17;
k = 0.
condition = SA[-1.5, 0.0, 0.0, -2.5, 0.0, 0.0]
p = SA[a, b, c, d,
        s, xr, r, I, xv, k1, k2, k]

pvalues = range(0, stop = 1, length = 500)
i = 4
plane = (1, 0.0)
tf = 100000.0
p_index = 12

output = produce_orbitdiagram(ds_HR, plane, i, p_index, pvalues,
                              tfinal = tf, Ttr = 50000.0;
                              diffeq = (alg = AutoVern9(Rodas5()), 
                                        abstol = 1e-11, reltol = 1e-11,
                                        dense = false,  maxiters = 10000000))

fig = Figure()
ax = Axis(fig[1,1]; xlabel = L"k", ylabel = L"x_2")
for (j, p) in enumerate(pvalues)
    scatter!(ax, fill(p, length(output[j])), output[j];
        color = ("black", 0.5), markersize = 0.7
    )
end
fig