There is an increasing literature on random non self-adjoint infinite matrices with motivations ranging from condensed matter physics to neural networks. Many of these operators fall into the class of ``pseudoergodic'' operators, which allows the elimination of probabilistic arguments when studying spectral properties. Parallel to this is the increased awareness that spectral properties of non self-adjoint operators, in particular stability, may be better captured via the notion of pseudospectra as opposed to spectra. Although it is well known that the finite section method applied to these matrices does not converge to the spectrum, it is often found in practice that the pseudospectrum behaves better with appropriate boundary conditions. We make this precise by giving a simple proof that the finite section method with periodic boundary conditions converges to the pseudospectrum of the full operator. Our results hold in any dimension (not just for banded bi-infinite matrices) and can be considered as a generalisation of the well known classical result for banded Laurent operators and their circulant approximations. Furthermore, we numerically demonstrate a convergent algorithm for the pseudospectrum including cases where periodic boundary conditions converge faster than the method of uneven sections. Finally we show that the result carries over to pseudoergodic operators acting on $l^p$ spaces for $p\in[1,\infty]$.