One thing to note is that, while you have a system of two equations, since each one is a second order equation, you have a fourth order system. That means that the phase space is 4D, so the separatrices are 3D manifolds in 4D space, and orbits are 1D manifolds (curves) in 4D space. You can project the orbits onto the plane, but they lose their property of being non-self-intersecting. Thus, most of the math that goes into this cannot be visualized, and so even the simplest coupled oscillators can provide a real challenge from a mathematical standpoint. Most of the theory on these things has been developed in the last 50 or so years, so keep that in mind when reading this. Also keep in mind that I'm a mathematician who studies dynamics, not a physicist or engineer, so while I have experience with numerical methods, my background is in analytical and theoretical methods.
For most strongly nonlinear systems, numerics are the only way (currently) to determine precisely the behavior, especially with respect to various parameters. For varying parameters, we usually refer to bifurcation theory. For a given set of parameters, we can identify the basins of attraction, which are the IC regimes where the system tends to each behavior. The boundaries between these regimes are called separatrices (singular: separatrix), or invariant manifolds. In autonomous systems, these manifolds are nice and simple (meaning not fractals) and each basin is a connected domain (two IC's in the same regime always have a curve connecting them entirely in the regime).
So, what we can do to reduce the trial and error to reasonable amounts is to try to find the fixed points and their linearized behavior, approximate the separatrices to get a partitioned space, then choose some sample IC's within each partition to try and determine the behavior of the system in that regime. You may stumble upon some points which generally move to fixed points, some which move to periodic orbits, some which move to quasiperiodic orbits, and some which move to chaos. Within the regime where a cycle happens, there must be some fixed point. You can transform that fixed point to the origin, then get the amplitude/phase system by assuming $y_i = r_i \cos(t+\psi_i)$ for $i=1,2$, then use variation of parameters and averaging to get an approximation of the cycle.
For the linear stability analysis, see any undergraduate book on Linear and Nonlinear ODE's. Approximation methods for invariant manifolds usually come in the form of some power series methods, but you can search for those. Once you've determined some approximate partitioning, just pick a bunch of initial conditions, then numerically solve the equations (using, say, ode45 with MATLAB) for each of the IC's and plot the results on the $(y_1,y_2)$-plane. For the variation of parameters and averaging, you can look at some books on Nonlinear Oscillations or Perturbation Methods, for example by Ali Nayfeh.
