The author presents an iterative algorithm to solve Toeplitz and non-Toeplitz block matrix equations. The development is based upon some well-known matrix iterative techniques. The algorithm is developed for the ideal case where the individual blocks in the autocorrelation matrix are Toeplitz, and it is then extended to a more general least squares data formulation case. It is shown that the algorithm requires fewer computations than the direct matrix inversion methods and is very simple to implement. The algorithm is applied to compute the spectral estimates of 2-D data of very small size based on the least squares data formulation with a quarter plane support.< >