A novel method of obtaining optimal reduced-order models of linear system transfer functions is presented. It uses the popular multipoint Padé approximation technique in an iterative way to generate efficiently the optimal models. Central to the method is a new way of calculating Padé approximants about many points by reducing them to equivalent Taylor series approximants. Optimal reduced-order models for impulse and step inputs are considered, and it is seen how the method may be extended to ramp and other polynomial inputs. Numerical examples are given to demonstrate the method.