Approximation of stable linear dynamical systems by means of so-called Laguerre and Kautz Functions, which are the Laplace transforms of a class of orthonormal exponentials, is studied. Since the impulse response of a stable finite dimensional linear dynamical system can be represented by a sum of exponentials (times polynomials), it seems reasonable to use basis functions of the same type. Assuming that the transfer function of a system is bounded and analytic outside a given disc, it is shown that Laguerre basis functions are optimal in a mini-max sense. This result is extended to the "two-parameter" Kautz functions which can have complex poles, while the poles of Laguerre functions are restricted to the real axis. By conformal mapping techniques the "two-parameter" Kautz approximation problem is recast as two Laguerre approximation problems. Thus, the well-developed theory of Laguerre functions can be applied to analyze Kautz approx approximations. Unilateral shifts are used to further develop the connection between Laguerre functions and Kautz functions. Results on H-2 and H-infinity, approximation using Kautz models are given. Furthermore, the weighted L(2) Kautz approximation problem is shown to be equivalent to solving a block Toeplitz matrix equation.