Your approach is putting the I and Q samples on a "sub-carrier" at frequency $f$. If that is really what you wanted to do, the typical demodulation approach is to multiply by the sub-carrier again and low pass filter to get I or Q. Dividing by cosine or sine creates time samples that go to infinity since you will divide by zero with those functions. The multiplication operation produces a baseband component and a component at twice the carrier frequency; this is made clearer with a quick refresh of the following trigonometric identities:
\begin{align}
\cos(\alpha)\cos(\beta) &= \frac{1}{2}\cos(\alpha-\beta)+\frac{1}{2}\cos(\alpha+\beta)\\
\sin(\alpha)\sin(\beta) &= \frac{1}{2}\cos(\alpha-\beta)-\frac{1}{2}\cos(\alpha+\beta)\\
\sin(\alpha)\cos(\beta) &= \frac{1}{2}\sin(\alpha-\beta)+\frac{1}{2}\sin(\alpha+\beta)
\end{align}
So in your case:
$$
\cos(2\pi ft)\cos(2\pi ft) = \frac{1}{2}\cos(0)+\frac{1}{2}\cos(2\pi ft) = \frac{1}{2} +\frac{1}{2}\cos(4\pi ft)
$$
The $\frac{1}{2}$ represents the amplitude of I and the $\frac{1}{2}\cos(4\pi ft)$ is the doubled frequency component that is removed with a low pass filter.
However I am not sure that you intended or need to have such a "sub-carrier"? If you were going to simulate QPSK modulation for example, you can modulate directly to your final carrier as in:
\begin{align}
I_t &= I \cos(2\pi f_c t), \quad\text{where}\quad I = {1,-1}\\
Q_t &= Q \sin(2\pi f_c t), \quad\text{where}\quad Q = {1,-1}
\end{align}
And further the I and Q samples are sent simultaneously by adding the modulated outputs $I_t$ and $Q_t$ above (hence the increased data rate; we send two bits in each "symbol).
To demodulate, assuming we have done the carrier recovery work to have a good estimate of $f_c$ with all phase offsets removed, we can multiply the received sequence by $2\cos(2\pi f_c t)$ and low pass filter to recover I, and multiply the received sequence by $2\sin(2\pi f_c t)$ and low pass filter to recover Q:
Noting the demodulation steps below where LPF() denotes removing the $2f_c$ frequency component with a low pass filter:
\begin{align}
I &= LPF\big\{2 I_t\cos(2\pi f_c t)\big\}\\
Q &= LPF\big\{2 Q_t\sin(2\pi f_c t)\big\}
\end{align}
The Q component goes to zero on the I output since the result is in the $\sin(\alpha)\cos(\beta)$ form resulting in $\sin(0) = 0$:
$$
0 = LPF\big\{2 Q_t\cos(2\pi f_c t)\big\}
$$
And similarly the I component goes to zero on the Q output:
$$
0 = LPF\big\{2 I_t\sin(2\pi f_c t)\big\}
$$
Therefore we can "multiplex" I and Q together given this ability to separate the two that operate simultaneously and independently in time.
Note I highly prefer to think and model using exponential frequency terms instead of sines and cosines - it is so much simpler! In fact this best explains implementations with I and Q paths as well, since to implement a complex signal you will need two real sequences $(e^{j\omega t} = \cos(\omega t) + j\sin(\omega t))$. Please see Frequency shifting of a quadrature mixed signal for more details on that.