Is the momentum operator well-defined in the basis of standing waves? 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?
 A: The Hilbert space for a particle in a box has no associated momentum operator (as a momentum operator implies that the state space is invariant under space translations). So your attempt do define one leads to artifacts.
[Edit] In general, the Hamiltonian must be a densely defined operator on the physical Hilbert space. But a square well potential is too singular as an operator on the space of all square integrable functions on $R$ (i.e., no densely defined subspace is mapped by it into the Hilbert space), hence its Hilbert space cannot be the standard Hilbert space. 
[Edit2] Note that the usual formula for the momentum is not a canonical momentum on the restricted space of wave functions defined only in the box, as it fails to satisfy either the commutation relation or does not preserve the boundary conditions. 
The appropriate replacement for the momentum operator is the mode counting operator $\hat n$ with $\hat n|n\rangle=n|n\rangle$; see also the derivation of Qmechanic.
A: Applying $\hat p$ to a Hamiltonian eigenfunction gives another state indeed, but this state can be decomposed in the eigenstates of the Hamiltonian. Numerically it is good everywhere except, maybe, the boundary points, so for integration purposes such a decomposition (superposition) is good. You did not obtain a diagonal matrix exclusively because you took a finite matrix for $\hat p$.
EDIT: Yes, at $t=0$ you can form a wave packet moving in a certain direction from "plane waves" $\exp(i k_n x)$, as usual, and then you can replace each momentum eigenfunction $\exp(i k_n x)$ with a linear combination of sines (Hamiltonian eigenfunctions $\psi_m (x)$ with the corresponding time exponentials $\exp(−i E_m t),\; t>0)$. In particular, if you choose the largest possible wave packet with a certain momentum $\psi(x,t=0)\propto \exp(-ik_n x),\;0<x<L$, then in later moments it will not keep its form because of reflections and the momentum will not conserve.
EDIT 2: Any standing wave $\psi(x)$ on an interval $(0,L)$ can be decomposed into a superposition of the momentum eigenfunctions $ \exp(\pm|k_n|x)$. In our case the superposition contains only two running waves with certain momenta: $\sin(k_n x)=(e^{ik_n x}-e^{-ik_n x})/(2i)$, so the momentum operator exists and there is no problem here.
A: A geometric approach can be achieve via an orbifold construction. 
I) Start with a free 1D quantum particle $\psi(x)=\psi(x+2\pi R)$ on a circle $S^1$ of circumference $2\pi R=2L$. The momentum $p=n\frac{\hbar }{R}$ is then quantized, $n\in\mathbb{Z}.$ The Hilbert space $H$ has momentum basis  $(| n\rangle)_{n\in\mathbb{Z}}$, where
$$\hat{p} | n\rangle~=~n\frac{\hbar }{R}| n\rangle, \qquad 
\langle n| m\rangle~=~\delta_{n,m}, \qquad n,m\in\mathbb{Z}. $$
II) The orbifold
$$I ~=~ S^1/ \mathbb{Z}_2~=~ S^1/\sim$$
is an interval of length $L$. A position $x \sim -x$ is identified with the opposite position $x\in S^1$ via a  discrete $\mathbb{Z}_2$-group action. To implement the two Dirichlet boundary conditions 
$$\psi(x=0)~=~0~=~\psi(x=L),$$ 
we impose that the wave function should be odd $\psi (x) = -\psi(-x)$. This kills half of the previous Hilbert space. The $\mathbb{Z}_2$-reduced Hilbert space $H/\sim~$ has basis  $(| n\rangle_{\sim})_{n\in\mathbb{N}}$, where
$$|n \rangle_{\sim}~:=~\frac{|n\rangle  - |{-}n\rangle}{\sqrt{2}}=-|{-}n \rangle_{\sim}, \qquad 
{}_{\sim}\langle n| m\rangle_{\sim}~=~\delta_{n,m}- \delta_{n,-m}, \qquad 
n,m\in\mathbb{Z}.$$
We can implement the orbifold construction through a parity operator $P$, 
$$P\hat{x}P~=~-\hat{x}, \qquad P\hat{p}P~=~-\hat{p}, \qquad P^2~=~{\bf 1}, \qquad P^{\dagger}~=~P.$$
A ket state $|\psi \rangle~\sim 0 $ is identified with zero if it is even, $P|\psi \rangle= |\psi \rangle$. The new basis $|n \rangle_{\sim}$ is odd, $P|n \rangle_{\sim}=-|n \rangle_{\sim}$. The absolute value $\hat{|p|}$ of the momentum operator is a $P$-even (unbounded) operator. It is diagonal in the $| n\rangle_{\sim}$ basis,
$$ \hat{|p|}~ | n\rangle_{\sim}~=~ |n|\frac{\hbar }{R}| n\rangle_{\sim},\qquad  {}_{\sim}\langle n| ~\hat{|p|}~|m\rangle_{\sim}~=~  |n| \frac{\hbar }{R}~ {}_{\sim}\langle n| m\rangle_{\sim}\qquad 
n,m\in\mathbb{Z}. $$
A: The operator is well-defined for instance in the position basis the operator corresponding to $\Psi(x) = \langle x | n \rangle$ is given by the ordinary form $\hat{p} = -i \hbar (d/dx)$ [*].
What happens is that $\Psi(x)$ are not eigenfunctions of the momentum operator $\hat{p}$ and, therefore, you cannot find a spectral decomposition of the operator in that basis (and, as a consequence, you cannot find a diagonal matrix).
Notice that $\Psi(x)$ are eigenfunctions of the operator $\hat{p}^2$ and it is possible to obtain the spectral decomposition
$$\hat{p}^2 = \sum p_n^2 | n \rangle  \langle n |$$
with the corresponding eigenvalues $p_n^2 = n^2h^2/4L^2$.
As a final note, it is not possible to obtain the matrix $\mathbf{p}$ from the matrix $\mathbf{p}^2$. First because the sign of $p$ is not defined and second because $\pm p$ does not correspond to a measurement of the momentum.
[*] By simplicity I am assuming one dimension.
