Compute the damping value for two masses in a harmonic oscillator Consider that this is a simple mass spring damper system:
$$ \ m \frac{d^{2}x}{dt} = F - b\frac{dx}{dt} - kx\ $$

What I allready know is the force $ F $ and the mass $ m $. Then I can find the spring stiffness k by sett all the derivatives to zero.
Then I can compute the damper $ b $ throught this formula:
$$ \ b = \zeta 2 \sqrt{k m}\ $$
According to Wikipedia, the damping ratio $ \zeta = 0.7 \ $ is a good number.
And this really works too! I have test it by my own ODE simulations, and it works OK! Perfect!
But the question is:
If I have two or more ODE:s with diffrent spring stiffness and dampness like this:
$$ \ m_{2}\frac{d^{2}z}{dt} = F - C_{2}(\frac{dz}{dt} - \frac{dy}{dt}) - k_{2}(z - y) \\ m_{1}\frac{d^{2}y}{dt} = C_{2}(\frac{dz}{dt} - \frac{dy}{dt}) + k_{2}(z - y) - C_1\frac{dy}{dt} - k_1 y\  $$
I can compute the stiffness $ k_1 , k_2 $ if I know the masses $ m_1 , m_2 $ and the force $ F $ , and also I need to set some position values for $ y , z $
Very easy! But how do I compute the damper $ C_1 , C_2 $ if I determine that $ \zeta = 0.7 $ ?

 A: 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}$.
A: In the general case there is no simple answer to this, because the vibration modes of the damped system are different from the undamped system. Even worse, the modes of the damped system may be complex - i.e. different degrees of freedom do not move in phase with each other.
But for your specific model, those problems don't apply. First write the equations in matrix form:
$$M \ddot x + C \dot x + Kx = F$$
where
$$\begin{align}
M &= \begin{bmatrix}m_1 & 0 \\ 0 & m_2 \end{bmatrix} \\[6pt]
C &= \begin{bmatrix}c_1 + c_2 & -c_2 \\ -c_2 & c_2 \end{bmatrix} \\[6pt]
K &= \begin{bmatrix}k_1 + k_2 & -k_2 \\ -k_2 & k_2 \end{bmatrix} \\[6pt]
x &= \begin{bmatrix}y \\ z \end{bmatrix}
\end{align}$$
For any values of the $m$'s, $c$'s, and $k$'s, we can express $C$ as a linear combination of $M$ and $K$, i.e. $C = \alpha M + \beta K$ for some constants $\alpha$ and $\beta$. This property of $C$ is known as Rayleigh damping and has the nice consequence that the mode shapes of the damped and undamped systems are identical.
So, you can solve your problem as follows:


*

*Find the modes of the undamped system $M\ddot x + Kx = F$

*Set the damping factors for each mode to your preferred value

*Find the damping matrix $C$

*Identify the physical coefficients $c_1$ and $c_2$.


There is a potential problem with this: you might find that $c_1 < 0$ which is not physically meaningful. But that is most likely to happen if you try to make the damping ratios for the two modes very different from each other, and apparently you don't want to do that.
Google will find plenty of references to Rayleigh damping, and it should be covered in a good dynamics textbook in the section on multi-degree-of-freedom (MDOF) systems.
In fact you don't need to go through all the steps I described above: from the general theory of Rayleigh damping, you can find the constants $\alpha$ and $\beta$ directly from the modal frequencies of the undamped system, and then get $c_1$ and $c_2$ directly. (But writing out all the theory behind that is more work than I'm prepared to do right now!)
