The hydrogen Hamiltonian is indeed precisely of the form I wrote. You just used the usual shorthand-notation not noting the identity operators for the single-body piece of the Hamiltonian. Of course, the entire system lives on the product space ##\mathcal{H}_{\text{p}} \otimes \mathcal{H}_{\text{e}}##.
It's however extremely unwise to use the naive product basis to solve the energy eigenvalue problem. As
@stevendaryl told you, it's much better to start by finding complete sets of compatible observables, including the Hamiltonian to find the basis where the diagonalization of ##\hat{H}## is easy. Since the observables in the complete set of compatible observables must commute with the Hamiltonian and among each others, you just have to look for the conserved quantities of your dynamical problem.
Here, it's indeed easy to find such a complete set of compatible observables by just thinking about the classical problem: It's a closed system of two Newtonian particles with a conservative central force. Thus the total momentum (center-of-mass momentum) is conserved, providing three compatible conserved observables (the three components of center-of-mass momentum).
The next step is to introduce new variables to express you Hamiltonian. You want the center-of-momentum components,
$$\vec{P}=\vec{p}_{\text{e}}+\vec{p}_{\text{p}}$$
Then it's good to introduce the relative position
$$\vec{r}=\vec{r}_{\text{e}}-\vec{r}_{\text{p}}$$
and it's canonical momentum ##\vec{p}##.
Expressing the Hamiltonian in these coordinates, you'll find after some simple algebra (using Heaviside-Lorentz units)
$$\hat{H}=\frac{\vec{P}^2}{2M} + \frac{\vec{p}^2}{2 \mu} -\frac{e^2}{4 \pi r}.$$
Here ##M=m_{\text{e}}+m_{\text{p}}## and ##\mu=m_{\text{e}} m_{\text{p}}/(m_{\text{e}} + m_{\text{p}})##.
Since ##\vec{P}## commutes with ##\vec{p}## and ##r## you can use ##P## eigenfunctions to diagonalize the Hamiltonian, i.e., you split your Hilbert space ##\mathcal{H}=\mathcal{H}_{\text{p}} \otimes \mathcal{H}_{\text{e}} = \mathcal{H}_{\text{CM}} \otimes \mathcal{H}_{\text{rel}}##.
Now you only need to solve for the effective one-body problem with the Hamiltonian
$$H_{\text{rel}}=\frac{\vec{p}^2}{2 \mu} -\frac{e^2}{4 \pi r}.$$
Since everything is radial symmetric, you can use ##H_{\text{rel}}##, ##\vec{\ell}^2##, and ##\ell_z## as a complete set of energy eigenvectors, where ##\vec{ell}=\vec{r} \times \vec{p}## is the orbital angular momentum of the relative motion.
So the appropriate basis to solve for the energy eigenvalue problem is ##|\vec{P},E_{\text{rel}},\ell,m \rangle##. The eigenspectra are ##\vec{P} \in \mathbb{R}^3##, ##\ell \in \{0,1,2,\ldots\}## and (for each given ##\ell##) ##m \in \{-\ell,-\ell+1,\ldots,\ell \}##. The Eigenvalues of ##\vec{\ell}^2## are ##\hbar^2 \ell(\ell+1)## and that of ##\ell_z## are ##\hbar m##.