For the specific case of the particle-in-a-box problem, the momentum operator is much trickier to handle than what you're allowing for, but that's not the issue at the heart of your current confusion.
If you know that two operators $A$ and $B$ commute, i.e. $[A,B]=0$, then there's two common ways to understand the consequences, one of which is right and one of which is wrong:
You are guaranteed the existence of at least one shared eigenbasis of both operators; but
You are not guaranteed that all eigenbases of $A$ will be eigenbases of $B$, or vice versa.
(For other occurrences of the misconception that the second point does happen, see e.g. this answer or this one.)
Thus, the 'contradiction' you observe is also present for the case of the free particle, without the box: the wavefunction $\psi(x)=\cos(kx)$ is an eigenfunction of $\hat H=\frac{1}{2m}\hat P^2$, but it is not an eigenfunction of $\hat P$. That is, the fact that the hamiltonian's eigenfunctions over a compact interval are not momentum eigenfunctions is not surprising and it is not specific to that configuration.
But, that said, you can still ask something like
well, OK, so some $\hat H$ eigenfunctions in the particle-in-a-box problem are not $\hat P$ eigenfunctions and that's not a problem, but $[\hat H,\hat P]=0$ is still true, so shouldn't I be guaranteed a shared eigenbasis, even if it isn't the one I started with?
and it's a reasonable question. Here what happens is that the subtleties with the compact-interval momentum operator kick in, at a much deeper level than just the commutation: the result, in full, reads
If two self-adjoint operators $A$ and $B$ commute, then there exists at least one shared eigenbasis,
and it breaks because $\hat P$-in-a-box is not a self-adjoint operator: it is symmetric, but it has domain problems that prevent it from being self-adjoint. The consequences of this are as deep as they get: there simply isn't a momentum eigenbasis in this Hilbert space. You can extend the momentum operator to make it self-adjoint; this extension is not unique, but there is a reasonable choice (setting $\alpha=0$ in V. Moretti's answer) that comes close to making $\hat H$ the square of the extended $\hat P_\alpha$, but ultimately this can't work, as they have different domains. (More specifically, the domain of $\hat H$ is contained in the domain of $\hat P_\alpha$, but the eigenfunctions of $\hat P_\alpha$ do not fall into that subspace.)
The double-slit experiment led us to think that the intensity of a classical wave is proportional to the probability of finding the particle there.
Consequently, the wavefunction is such that its squared modulus represents the probability
$||\psi(x_0) || ^2 dx $ represents the probability of finding the particle at $x=x_0$ (in an infinitesimal environment of $x_0$)
The integral of that quantity
$$ \int_a^b ||\psi(x) || ^2 dx $$
represents the probability of finding the particle in the interval $a,b$. It is usually written as
$$ \int_a^b \psi^*(x)\cdot\psi(x)\cdot dx $$
which is the same.
The integral along the whole space means, obviously, the probability of finding it anywhere, so it must be 1, that's why $\psi$ mus be "normalized".
Finally, the $1s$ orbital is just the solution of the equation. It turns out to be like that. You have to solve an ideal hydrogen atom, that is, a Coulomb's potential. You find out that the radial functions of hydrogen are like that.
However, probability being maxium does not mean that the electron is always there.
Best Answer
In your example the difficulty is that you're taking linear combination of eigenstates of $p$ but with different eigenvalues, so the resulting combination is no longer an eigenstate of $p$, even if the pieces are separately eigenstates.
An alternate example would be the simple case $[\hat H,\hat L^2]=0$ and a hydrogen atom state with $n=2$ so that $\ell=0,1$ can occur. Then $\{\vert n\ell m\rangle\}$ are simultaneous eigenvectors of $\hat H$ and $\hat L^2$ but a combination of these containing different $\ell$s will not be a simultaneous eigenstate of both since different $\ell$ states have different eigenvalues of $\hat L^2$.