Rather than the frequency domain, let's look at this in the time domain and particularly, the characteristic equation associated with a linear homogeneous 2nd order differential equation for some system:
\$r^2 + 2 \zeta \omega_n r + \omega^2_n = 0\$.
If the roots of the characteristic equation are real (which is the case if \$\zeta \ge 1\$), the general solution is the sum of real exponentials:
\$Ae^{\sigma_1 t} + Be^{\sigma_2t} \$
where
\$\sigma_1 = -\zeta \omega_n + \sqrt{(\zeta ^2 - 1)\omega^2_n} \$
\$\sigma_2 = -\zeta \omega_n - \sqrt{(\zeta ^2 - 1)\omega^2_n} \$
Since these are real exponentials, there is no oscillation in these solutions.
If the roots are complex conjugates (which is the case if \$\zeta < 1\$), the general solution is the sum of complex exponentials:
\$e^{\sigma t}(Ae^{j\omega t} + Be^{-j\omega t})\$
where
\$\sigma = -\zeta \omega_n\$
\$\omega = \sqrt{(1 - \zeta ^2)\omega^2_n}\$
This solution is a sinusoid with angular frequency \$\omega\$ multiplied by a real exponential. We say the system has a "natural frequency" of \$\omega\$ for a reason that I think is obvious.
Finally, setting \$\zeta = 0\$ (an undamped system) , this solution becomes:
\$Ae^{j\omega_n t} + Be^{-j\omega_n t}\$
which is just a sinusoid of angular frequency \$\omega_n\$.
In summary, a system may or may not have an associated natural frequency. Only systems with \$\zeta < 1\$ have a natural frequency \$\omega\$ and only in the case that \$\zeta = 0\$ will the natural frequency \$\omega = \omega_n\$, the undamped natural frequency.
You want to use the following property of Laplace transform:
$$\mathscr{L}\left({\frac{dx(t)}{dt}}\right)(s) = s\mathscr{L}\left(x(t)\right)-x(0)$$
This allows you to easily move between differential equations and polynomial equations.
Time domain to frequency domain: Take your first equation for example
$$\frac{d^2 x(t)}{dt^2} + 2 \zeta \omega_n \frac{dx(t)}{dt} + \omega_n^2 x(t) = f(t).$$
If we denote the Laplace transform of \$x(t)\$ by \$X(s)\$ and \$f(t)\$ by \$F(s)\$ and apply Laplace transform to this equation then this property implies
$$s^2 X(s) + 2 \zeta \omega_n sX(s) + \omega_n^2 X(s) = F(s)$$
where for simplicity I assume that \$x(0) = f(0) = 0\$.
The transfer function is defined as the ratio:
$$\frac{X(s)}{F(s)} = \frac{1}{s^2+2 \zeta \omega_n s + \omega_n^2}$$
Frequency domain to time domain: Lets try the example \$G(s)=\frac{1}{s^3+1}\$ then we have by definition that
$$G(s)F(s) = X(s)$$
which implies
$$F(s) = s^3X(s)+ X(s).$$
Taking inverse Laplace transform we find that
$$f(t) = \frac{d^3 x(t)}{dt^3} + x(t).$$
Hopefully this allows you to see the pattern in general.
Best Answer
I think you mainly need help understanding what it means for the system to be ‘second-order dominated on the verge of instability.’
It means that there is a pair of complex poles just barely to the left of the imaginary axis (or perhaps on the axis, which is marginally stable), with no other poles that are also close to the imaginary axis.