Settling time is usually defined as when when the step response is within 2% of the steady state response. Even though the difference between the step response and the steady state response gets bigger for all times when you use a bigger gain $K$, the percentage difference remains the same. That is why $K$ does not affect settling time.
The poles of your system are,
$$
p = \omega_n \left(-\zeta \pm \sqrt{\zeta^2 - 1}\right).
$$
The transient response of such a system will look like,
$$
y_t(t) = C_1\, e^{\omega_n \left(-\zeta + \sqrt{\zeta^2 - 1}\right) t} + C_2\, e^{\omega_n \left(-\zeta - \sqrt{\zeta^2 - 1}\right) t},
$$
where $C_1$ and $C_2$ depend on the initial conditions.
For $\zeta<1$ the poles will be each others complex conjugate and the transient can also be written as,
$$
y_t(t) = e^{-\omega_n\, \zeta\, t} \left(C_3 \cos\left(\omega_n\,\sqrt{1 - \zeta^2}\, t\right) + C_4 \sin\left(\omega_n\,\sqrt{1 - \zeta^2}\, t\right)\right).
$$
For a step response the total response can be written as the sum of the steady state and the transient response ($y(t)=y_{ss}(t)+y_t(t)$), with the steady state equal to $K$ (which can be found by setting $s$ to zero). For a step response it is assumed that the system is initial at rest, so position and velocity are zero for all $t\leq0$. Applying these constraints yields $C_3=-K$ and $C_4=0$,
$$
y(t) = K \left(1 - e^{-\omega_n\, \zeta\, t} \cos\left(\omega_n\,\sqrt{1 - \zeta^2}\, t\right)\right).
$$
Using an upper bound of one for the cosine, then the settling time can be approximated by only considering the exponential,
$$
\left|\frac{y(t)-y_{ss}(t)}{y_{ss}(t)}\right| \approx e^{-\omega_n\, \zeta\, t}.
$$
It can be noted that if you work out the left hand side of this equation yourself you will see that all $K$'s cancel.
The settling time for this approximation can be found by solving when $e^{-\omega_n\, \zeta\, t}=0.02$. This has the solution,
$$
T_{settling} = \frac{-\ln(0.02)}{\zeta\,\omega_n},
$$
$-\ln(0.02)$ is indeed close to $4$.
Solving for the constants for two real poles (when $\zeta\geq 1$) you get,
$$
C_1 = K \frac{\sqrt{\zeta^2 - 1} + \zeta}{2 \sqrt{\zeta^2 - 1}},
$$
$$
C_2 = K \frac{\sqrt{\zeta^2 - 1} - \zeta}{2 \sqrt{\zeta^2 - 1}}.
$$
For $\zeta$ equal or slightly bigger than 1 you will have to solve the following equation for $t$ in order to get the settling time,
$$
\left|e^{-\omega_n\, \zeta\, t} \left(\cosh\left(\omega_n \sqrt{\zeta^2 - 1}\, t\right) + \frac{\zeta}{\sqrt{\zeta^2 - 1}} \sinh\left(\omega_n \sqrt{\zeta^2 - 1}\, t\right)\right)\right| = 0.02.
$$
When $\zeta\gg1$ then $\zeta\, (\zeta^2 - 1)^{-1/2} \approx 1$ and $\cosh(t)+\sinh(t)=e^t$ can be used. Solving the above equation then has the solution,
$$
T_{settling} = \frac{-\ln(0.02)}{\left(\zeta - \sqrt{\zeta^2 - 1}\right) \omega_n}.
$$
This can also be seen as taking the biggest (least negative) pole, instead of the real part of the complex conjugate poles. Because when $\zeta$ is large the difference between the two real poles is large and only the "slowest" pole will dominate the transient response for the settling time.
Your approach seems correct. Namely, near a breakaway point two real poles should "meetup" and become a pair of complex conjugate poles or vice versa, thus at a breakaway point one would have repeated (real) poles. Repeated poles also means that the derivative of the characteristic equation with respect to $s$ at the value of the pole should be equal to zero. For example define the general openloop as
$$
H_{ol}(s) = k \frac{N(s)}{D(s)}, \tag{1}
$$
with $N(s)$ and $D(s)$ polynomials in $s$. Thus the characteristic equation would be
$$
D(s) + k\,N(s) = 0. \tag{2}
$$
The derivative of the characteristic equation would be
$$
D'(s) + k\,N'(s) = 0, \tag{3}
$$
with $X'(s)$ denoting the derivative of $X(s)$ with respect to $s$. However, inorder to assure that $s$ is a pole of $(2)$ one can indeed use
$$
k = -\frac{D(s)}{N(s)}. \tag{4}
$$
It can be noted that the $k$ in $(2)$ should not change value and thus do not contribute additional terms to the derivative in $(3)$. Substituting $(4)$ into $(3)$ yields
$$
D'(s) - \frac{D(s)}{N(s)}N'(s) = \frac{D'(s)\,N(s) - D(s)\,N'(s)}{N(s)} = 0. \tag{5}
$$
Since $N(s)$ is a polynomial $(5)$ is equivalent to
$$
D'(s)\,N(s) - D(s)\,N'(s) = 0. \tag{6}
$$
Even though it is stated that $k$ should remain constant, it can be shown that $(6)$ is equivalent to setting the derivative of $k$ with respect to $s$ to zero. Namely, by using the quotient rule it can be shown that the derivative of $k$ is equal to
$$
k' = -\frac{D'(s)\,N(s) - D(s)\,N'(s)}{N(s)^2}. \tag{7}
$$
When using again that $N(s)$ is a polynomial yields that setting $(7)$ to zero is equivalent to $(6)$.
Substituting your openloop transfer function does indeed yield
$$
2s^3+18s^2+48s+38=0. \tag{8}
$$
The only thing I could remark on is that you incorrectly rounded the corresponding roots of $(8)$. For example the first root rounded to two decimal places should be $s=-1.47$.
It can be noted that the other roots of $(8)$ are also breakaway points. Only, if you substitute those values into $(4)$ yields negative values for $k$, while normally a rootlocus plot considers only $k\in [0,\infty)$.
Best Answer
Intuitively you can think that open loop poles tend to "push" and zeros tend to "pull" the root loci. Suppose you have only two real poles, then roots break away at exactly the halfway between them as in System 2 in your link. If you add another real pole as in System 3, it pushes the roots, hence the break away point.
Suppose you have real poles at -1, -3 and a zero at -5. Then zero tends to pull the roots and in fact a circle with center at -5 and radius a little less than 3 will occur. Hence, the break away point will be between -1 and -3, but closer to -3 due to the effect of the zero. Of course there will be another break away point close to -8.
In System 4, zero at -4 is closer to other poles than the pole at -5, hence it pulls the roots "more strong" then the pole pushes. But their effects will even out at infinity, hence asymptote has a 90 degree angle. Note that angle of the asymptotes are directly determined by the difference between number of roots and zeros.