In this paper, we propose a new model-order reduction technique for linear dynamic systems. The idea behind this technique is to transform the dynamic system function from the s-domain into the z-domain via the bilinear transformation, then use Prony?s [1, 2] or Shank?s [3] least-squares approximation methods instead of the commonly employed Pad? approximation method, and finally transform the reduced system back into the s-domain using the inverse bi-linear transformation. Simulation results for large practical systems show that this technique based on Prony?s and Shank?s methods give much higher accuracy than the traditional Pad? method, and result in lower-order approximations with negligible increase in simulation time.