You've designed a second-order band pass filter, so the number of terms in the denominator is just fine (it's a second-order polynomial). Your design is OK, apart from the scaling, i.e., the maximum of the magnitude of the transfer function is not equal to $1$, which might be desirable.
Your transfer function has the form
$$H(z)=\frac{z^2-1}{(z-p)(z-p^*)}=\frac{z^2-1}{z^2-2zr\cos(\omega_p)+r^2}\tag{1}$$
where $p$ is a complex-valued pole:
$$p=re^{j\omega_p}\tag{2}$$
with radius $r$ and pole angle $\omega_p$. Note that in general the pole angle $\omega_p$ and the peak frequency $\omega_0$ (where the magnitude of the frequency response $H(e^{j\omega})$ achieves its maximum) are not equal. However, in your case they are almost identical because $\omega_p$ is relatively close to $\pi/2$. See this answer for more details on this matter.
If you want the frequency response peak value to be equal to $1$, you need to scale your transfer function by the factor $(1-r^2)/2$ (also explained in the answer quoted above).
The pole angle $\omega_p$ depends on the desired center frequency, so it can't always be equal to $\pi/2$ as you seem to have thought (even though in your calculations you chose it according to the desired center frequency).
The Matlab function freqz computes the complex frequency response of a discrete-time filter.
For plotting the poles and zeros you can use the function zplane.
