Suppose I want to describe an arbitrary state of a quantum particle in a box of side $L$. The relevant eigenmodes are those of standing waves, namely $$ \left<x|n\right>=\sqrt{\frac{2}{L}}\cdot \sin \left(\frac{n\pi x}{L}\right)$$ In this basis, the operate $\hat{p}^2$ is diagonal by construction, so every state has a definite energy.
But suppose I want to construct this operator from the matrix elements of $\hat{p}$? These are ($\hbar=1$) $$\left<m|\hat{p}|n \right>=-\frac{2imn((-1)^{m+n}-1)}{L\cdot(m^2-n^2)}$$ which vanishes for $m$ and $n$ with the same evenness, making this an off-diagonal "checkered" matrix.
I tried taking a finite portion of this matrix and squaring it in order to obtain $\hat{p}^2$, but from what I could tell, this yields an on-diagonal checkered matrix, contrary to the construction that it would be strictly diagonal.
I assume this has something to do with the fact that the wavefunctions obtained by applying $\hat{p}$ do not satisfy the same boundary conditions as the basis functions (do not vanish). Would this problem be solved I had taken the the infinite matrix in its entirety before squaring? Or is this entire problem (decomposing $\hat{p}$ in this basis) ill-defined?