Hi Z^2,
Now that you have the K, C and M matrices, you can create a matrix equation to find the natural resonant frequencies. Generalizing to n masses instead of 3, Let
xv = [x;v]
be a 2nx1 column vector of n displacements and n velocities; and let the system have an overall time dependence of exp((g+i*w)*t). Then a time derivative gives a multiplicative factor of (g+i*w) and
[v;a] = d/dt [x;v] = [x;v]*(g+i*w)
For the eqns of motion you can arrive at the 2nx2n matrix equation
[I 0;0 M]*[x;v]*(g+i*w) = [0 I;-K -C]*[x;v]
Here Mfull = [I 0;0 M] is a block matrix of four nxn matrices, I is the nxn identity, and 0 is an nxn block of zeros, and similarly for [0 I;-K -C]. Ordinarily you would have to multiply both sides by the inverse of Mfull to obtain an eigenvalue equation. But Matlab provides an option in eig so that you don't have to actually use the inverse of Mfull:
Mfull = [I 0;0 M];
A = [0 I;-K -C]
[xv lambda] = eig(A,Mfull)
The lambdas are the natural frequencies of the system, lambda = g+i*w. The values of g had better come out negative, meaning exponential decay rather than exponential growth.
I assumed here that you were not applying any external forces at a specific frequency and looking for a steady state solution. If you were, then the solution is different, and easier.
