Convert low pass continuous time filter design to bandpass, discrete time
There is an error in the analog bandpass filter transfer function. If you apply the LP-BP transformation as you stated it you should get the following transfer function:

$$H(s)=\frac{s^2(\Omega_2-\Omega_1)^2}{s^4+\sqrt{2}(\Omega_2-\Omega_1)s^3+(\Omega_1^2+\Omega_2^2)s^2+\sqrt{2}\Omega_1\Omega_2(\Omega_2-\Omega_1)s+\Omega_1^2\Omega_2^2}$$

With the given band edge frequencies (note that the lower band edge of the discrete-time filter must be $2000\,\text{Hz}$, not $4000\,\text{Hz}$) you correctly got $\Omega_1=2$ and $\Omega_2=4$. Applying the bilinear transform to the analog bandpass transfer function with these band edges finally gives the following discrete-time transfer function:

$$H(z)=\frac{z^4-2z^2+1}{(10+3\sqrt{2})z^4+(12+2\sqrt{2})z^3+20z^2+(12-2\sqrt{2})z+10-3\sqrt{2}}$$