If eig has been working well for the size of your problem you could consider using
[U, D] = eig(K, M); %this solves K*U = M*U*D and would match eig(M\K) if M is invertible.
instead of computing the inverse of M. This will work safely also for singular matrix M. You can of course flip the order of passing in K and M, this only affects if the eigenvalues computed need to be inverted (is it K*x = lambda*M*x or M*x = lambda*K*x).
By convention, usually people will use the mass matrix as the second matrix if it makes sense to call one of the input a mass matrix.
For eigs, you can use the same syntax eigs(K, M) but here it's much more relevant than with EIG which eigenvalues you want (eigs computes only a subset, typically closest to zero or furthest away from zero). Depending on this, either K or M still need to be inverted inside of eigs.
eigs has quite a lot of special-case treatment for the case when the second matrix is symmetric positive definite, so if one of your matrices is that, it may be a good idea to pass it in as the second argument.
