To represent an harmonically vibrating continuous system having an infinite number of degrees of freedom by means of mathematical models of a finite number of co-ordinates, the relation between the excitation force vector and response displacement vector (the dynamic stiffness matrix [ D(λ)]) is inevitably frequency dependent. It is often required to compute the dynamic flexibility matrix [ Z(λ)], which is the inversion of [ D(λ)]. A new method is presented here to express [ Z(λ)] in terms of the eigensolutions of [ D(λ)]. It is similar to the classical format except that [ D(λ)] is no longer required to be expressible as a matrix polynomial in λ with constant coefficient matrices. In contrast to the state variables, which tend to increase the order of the matrices, only physical co-ordinates are concerned. [ D(λ)] may be non-symmetric and defective.