Finding the frequency response of a filter defined as a Z-domain transfer function
The analytic way is to substitute the variable $z$ by $e^{j\omega}$ to get the frequency response $H(\omega)$ (with $\omega = \frac{2 \pi f}{F_s}$) - that is to say, the frequency response is the $z$ transform evaluated on the unit circle.

Note that matlab has a built-in function for plotting the frequency response straight from filter coefficients (freqz), whose source code can be seen. Have a look! The implementation does not even involve solving the difference equation or filtering an impulse - you can simply divide the FFT of the numerator sequence by the FFT of the denominator sequence - padded up to the desired number of frequency points.