The purpose of this paper is to investigate methods for the iterative computation of approximately balanced reduced order systems. The resulting approach is both completely automatic once an error tolerance is specified and yields an error bound. This is to be contrasted with existing projection methods, namely PVL (Padé via Lanczos) and rational Krylov, which do not satisfy these properties. Our approach is based on the computation and approximation of the cross gramian of the system. The cross gramian is the solution of a Sylvester equation and therefore some effort is dedicated to the study of this equation with some new insights. Our method produces a rank $k$ approximation to this gramian in factored form and thus directly provides a reduced order model and a reduced basis for the orignal system. It is well suited to large scale problems because there are no matrix factorizations of the large (sparse) system matrix. Only matrix-vector products are required.