A linear system with a Toeplitz matrix is solved by the preconditioned conjugate gradient squared (PCGS) iterative method. A product of a band Toeplitz and a circulant matrix is used as preconditioner. Using a careful numerical approximation scheme to the generating function of the Toeplitz matrix, superlinear convergence is achieved, even in cases when the generating function is complex-valued and has multiple zeros. The method is applied to solve stationary probability distributions of Markovian queueing networks with batch arrivals, as well as to specially devised numerical test cases.