You are right to be suspicious of an answer that contradicts your intuition - but in this case the solution is right
The key (not clearly explained) is that you apply a constant force that is greater than the force needed to move one block: consequently that block will accelerate. As the spring compresses, the first block will accelerate less, and eventually it will decelerate. When it stops, the compression of the spring is maximum - and it applies a greater force to both blocks than the force you used to compress it.
So the solution is correct. If you still have difficulty believing this, think about a hammer: when you swing it, you apply a little bit of force over a long distance; when it hits the nail, it feels a lot of deceleration over a short distance. The force applied to the nail is greater than the force you applied to the hammer. Similar (not same) concept... You do work over a large distance with a little bit of force, to allow you to do the same work over a short distance (with a large force)
Does that help?
I am going to tackle a special case of the general problem stated in the question. I hope this would leads to some insight into how to handle the general problem.
Consider the 2-DOF system, but with each mass equal to $m_1=m_2 =\frac{m}{2}$ (so the total mass is $m$), each stiffness $k_1=k_2 = 2 k$ (so the combined stiffness is $k$) and each damping is $c_1=c_2 =2 c$ for the same reason.
To solve the problem, I formulated the equations of motion in terms of the absolute displacements $x_1$ and $x_2$, but then I switched the DOF of the system to relative displacements $ q_1 = x_1 $ and $q_2 = x_2-x_1$.
As a system of equations this is
$$ \pmatrix{\ddot{q}_1 \\ \ddot{q}_2 } + \left[ \matrix{ \frac{4 k}{m} & -\frac{4 k}{m} \\ -\frac{4 k}{m} & \frac{8 k}{m} } \right] \pmatrix{q_1 \\ q_2} + \left[ \matrix{ \frac{4 c}{m} & -\frac{4 c}{m} \\ -\frac{4 c}{m} & \frac{8 c}{m} } \right] \pmatrix{\dot{q}_1 \\ \dot{q}_2}= \pmatrix{0\\0} $$
This system has an "in-phase" solution of
$$\pmatrix{q_1 \\ q_2 } = \pmatrix{Q_1 \\ Q_2} \exp(-\beta t) \sin(\omega t) $$
Since this is a 2-DOF system it has two natural frequencies
$$ \omega_n^2 = \begin{cases}
\frac{2 k( 3-\sqrt{5})}{m} & \mbox{mode A} \\
\frac{2 k( 3+\sqrt{5})}{m} & \mbox{mode B} \end{cases} $$
These are found from the eigenvalues of the 2×2 matrix multiplying $\pmatrix{q_1 & q_2}$ in the equations of motion.
The values of $\beta$ and $\omega$ that solve the system are
$$\begin{aligned}
\beta & = \frac{c\, \omega_n^2}{2 k} \\
\omega^2 & = \omega_n^2 \left(1- \frac{c^2\, \omega_n^2}{4 k^2} \right)
\end{aligned} $$
Using the values of $\beta$ and $\omega$ that solve the system of equations, one can parametrise them as $$ \begin{aligned} \beta & = \omega_n \zeta \\ \omega & = \omega_n \sqrt{1-\zeta^2} \end{aligned}$$ which is exactly what you do in the 1-DOF system. This means that the stiffness $k$ and damping $c$ needed for a response with natural frequency $\omega_n$ and damping ratio $\zeta$ are
$$\boxed{ \begin{aligned}
k & = \frac{m \, \omega_n^2}{2 (3 \pm \sqrt{5})} \\
c & = \frac{m \, \omega_n \zeta}{(3 \pm \sqrt{5})}
\end{aligned} }$$
Also note that to target a specific half-life of $t_H$ where $\exp(-\beta t_H) = \frac{1}{2}$ you need a damping ratio of $\zeta = \frac{ \ln(2)}{\omega_n\, t_H}$.
Finally, for completeness I am going to mention that the response amplitudes are arbitrary, but must be related to each with $Q_2 = -Q_1 \frac{1\pm\sqrt{5}}{2}$.
Best Answer
Yes, one can calculate an equivalent mass. Take the following coordinates, that describe the position of the particles' centers with respect to the center of mass of the set, $G$. In the sketch that follows I have arbitrarily chosen $m_1 < m_2$, and so $G$ falls closer to $m_2$:
Since $G$ will remain fixed in the absence of external forces, it makes sense to use the coordinates above. The equations of motion of the two particles in these coordinates read $$m_1\frac{\mathrm{d}^2 x_1}{\mathrm{d} t^2} + c \frac{\mathrm{d}l}{\mathrm{d} t} + kl = 0 \\ m_2\frac{\mathrm{d}^2 x_2}{\mathrm{d} t^2} + c \frac{\mathrm{d}l}{\mathrm{d} t} + kl = 0$$ where $l = x_1 + x_2$ (the distance between centers). Adding these two equations and deviding through by $2$ one obtains $$\frac{m_1 + m_2}{2}\frac{\mathrm{d}^2 l}{\mathrm{d} t^2} + c \frac{\mathrm{d}l}{\mathrm{d} t} + kl = 0$$ which is of the form of the original equation (the one from Wikipedia) if one takes $$m = \frac{m_1 + m_2}{2}$$
Thus, the equivalent mass is just the arithmetic average of the two masses.