Consider a finite dimensional linear system perturbed by nonlinear perturbations, e.g. in 3D, consider a system of the form $$\dot{x}=f_x(x,y,z)+\varepsilon g_x(x,y,z),$$ $$\dot{y}=f_y(x,y,z)+\varepsilon g_y(x,y,z),$$ $$\dot{z}=f_z(x,y,z)+\varepsilon g_z(x,y,z)$$ where $f_x,f_y,f_z$ are linear functions, and $g_x,g_y,g_z$ are nonlinear functions, of $x,y,z$. We know that the unperturbed system cannot be chaotic, and that the perturbed system certainly can be chaotic. The question is how can I place some system specific arguments to determine whether it is chaotic or not, based on, for example, the value of the perturbation parameter $\varepsilon$ and the functional form of the nonlinear components. For example, consider $$\dot{x}=-x+z+\varepsilon xy,$$ $$\dot{y}=-y+\varepsilon zy^3,$$ $$\dot{z}=-z+\varepsilon xy^3.$$ Is it possible to determine, at least in terms of a lower bound, for what value of $|\varepsilon|\ll1$ this system becomes chaotic, if at all?
Is it certainly chaotic for $\varepsilon\equiv 1$?