Modal Parameters from dynamic system (numerical errors)

I have a dynamic system modeled by the Mass, Damping and Stiffness matrices [M], [C] and [K], respectively, i.e.,

[M]{x''} + [C]{x'} + [K]{x} = {0} (no exciting force, for natural vibration).

I want to estimate the modal parameters natural frequency wn and damping factor z for this system. I can't use the damp function because this system has an overdamped mode (z>1), which means damp returns a wrong answer (it states simple z=1, which is not true).

The way I know to get these parameters is to use the state space model, i.e.,

A = [zeros(n,n) eye(n,n)
     -inv(M)*K -inv(M)*C];
l = eig(A)

I know from theory that each pair of eigenvalues is of the form -wn(z+-sqrt(z^2-1)), which means I can write a code for each case (z is greater than or less than one), so I can obtain my modal parameters from the eigenvalues l=eig(A).

The problem is: the matrices M, C and K are obtained from a finite element discretization, which means that, for large matrices M, C and K, the inverse calculations aren't very good and A has a condition number to high, so the values of wn varies a lot with the size of the matrices. Is there any way to calculate these parameters without numerical errors?

PS: I can calculate the natural frequencies from

eig(K,M)

which gives me a (very) good answer. However, eig(K,M) does not give me the damping factor. Is there a way to calculate the damping factor? That is the parameter I am looking for, since eig(K,M) already gives me a good answer.

PS2: if the approach I'm giving to the problem is good (calculating the state space matrix A), then probably the problem is the numerical erros when calculating -inv(M)*K and -inv(M)*C, for large matrices M, C and K. Is there a way to do this calculations with minor numerical errors?