How do I find the phase response of a filter from the filter equation?
You need to determine the frequency response of the filter. The phase of this complex function is the filter's phase response. If you have the filter's difference equation, you need to transform it to the $\mathcal{Z}$-transform domain, where a delay of $k$ (with integer $k$) corresponds to a multiplication by $z^{-k}$:

$$Y(z)=c_1X(z)+c_2Y(z)z^{-1}\tag{1}$$

From $(1)$ you can determine the transfer function

$$H(z)=\frac{Y(z)}{X(z)}=\frac{c_1}{1-c_2z^{-1}}\tag{2}$$

The frequency response is the transfer function evaluated at $z=e^{j\omega}$, where $\omega$ is the normalized frequency in radians:

$$\omega=\frac{2\pi f}{f_s}\tag{3}$$

with $f_s$ the sampling frequency, and $f$ the actual frequency in Hz.

From $(2)$, the filter's frequency response is

$$H(e^{j\omega})=\frac{c_1}{1-c_2e^{-j\omega}}=|H(e^{j\omega})|e^{j\phi(\omega)}\tag{3}$$

where the phase response $\phi(\omega)$ is given by

$$\phi(\omega)=\text{atan2}\left(\frac{-c_2\sin(\omega)}{1-c_2\cos(\omega)}\right)\tag{4}$$

For more complicated filters it's usually easier to determine the frequency response (and the phase response) numerically by using an FFT. The (sampled) frequency response is obtained by element-wise division of the FFT of the numerator coefficients of the transfer function by the FFT of the denominator coefficients. The (sampled) phase is then obtained as the argument of this complex-valued vector. In Matlab/Octave, the frequency response can be obtained directly from the filter coefficients using the command freqz.