Understanding where does the second (stochastic) attractor of the system come from.

Modern ODE solver packages have events and actions. Events test for some condition, actions change state, parameters or event the ODE function at such points. The scipy solver solve_ivp only implements events and the hard-coded actions to register the event points (always) and to terminate the integration (by a flag).

ODE function and event can be implemented as

def model(t,u,a):
    x,y,z = u
    dx = x*(a - 0.06*x - y/(x+10))
    dy = y*(-1+2*x/(x+10)-0.405*z/(y+10))
    dz = z*z*(0.038-1/(y+20))
    return dx,dy,dz

a=1.803

u0 = [29.9,0.2,15.0]

def event(t,u,a): return u[1]-24
event.direction = -1

T = 40
res = solve_ivp(model,(0,10*T),u0, args=(a,), events=event, atol=1e-6, rtol=1e-9)
res = solve_ivp(model,(0,20*T),res.y_events[0][-1], args=(a,), events=event, atol=1e-6, rtol=1e-9)

print(res.t_events[0])

As we are looking for stable attractors, it makes sense to let the solver run for some time to get close to an attractor. The period is about $T=40$, letting it run for 20 periods may seem excessive, 10 are also sufficient.

Now loop through some random initial points from $[25,30]×[0,10]×[0,20]$ and plot the stabilized trajectories, and the event locations in the x-z-plane. One finds that the events lie on a straight line, thus also draw the path of the $x$ coordinate in a successor map.

fig = plt.figure(figsize=(7,4.5), constrained_layout=True)
gs = fig.add_gridspec(2, 3)
ax = [ fig.add_subplot(gs[k, 2]) for k in [0,1]]
ax3 = fig.add_subplot(gs[:, :2])
for k in range(10):
    u0 = np.random.uniform(size=3)*np.array([5,10,20])+np.array([25,0,0])
    res = solve_ivp(model,(0,10*T),u0, args=(a,), events=event, atol=1e-6, rtol=1e-9)
    res = solve_ivp(model,(0,10*T),res.y_events[0][-1], args=(a,), events=event, atol=1e-6, rtol=1e-9)
    ax3.plot(res.y[0],res.y[2], lw=0.8)
    ax[0].plot(res.y_events[0][:,0],res.y_events[0][:,2],'-+', lw=0.5)
    ax[1].plot(res.y_events[0][:-1,0],res.y_events[0][1:,0],'-+', lw=0.5)
for ak in [ax3,*ax]: ak.grid(); 
ax3.set_xlabel("$x$"); ax3.set_ylabel("$z$")
ax[0].set_xlabel("$x_E$"); ax[0].set_ylabel("$z_E$")
ax[1].set_xlabel("$x_E[k]$"); ax[1].set_ylabel("$x_E[k+1]$")
plt.show()

In the resulting plot the stable loop is well-isolated, the point in the lower left corners, and the alternating connection of the two loops of the chaotic attractor is also visible in the second and last plot.

enter image description here