This is not a solution, but some hopefully useful notes that may lead to other things useful towards a genuine solution, and do not fit into a comment.
I mentioned in a comment that the model should be exactly solvable. The issue is that your coefficients A and C are time-dependent and things get complicated. If there is any way to approximate/replace them by time-independent values, everything simplifies a lot, to the point where the problem becomes algebraically tractable. Simple oscillatory exponentials might also work to some degree.
In any case, the general idea for an exact solution is that the master equation's generator is of a form ${\mathcal L} = {\mathcal L}_0 + {\mathcal L}_1$, where ${\mathcal L}_0$ and ${\mathcal L}_1$ yield relatively simple nested commutation relations (as superoperators), which in turn allow a suitable disentangling of the time evolution (for either ${\mathcal L}$ or ${\mathcal L}^\dagger$) into simpler tractable factors. The time evolution of the number operator and its average may then be calculated exactly. An approach developed along these lines for the $C=0$ case and based on a su(1,1) algebra is given in http://lanl.arxiv.org/abs/0710.2724. Details on their su(1,1) method in http://lanl.arxiv.org/abs/quant-ph/0112090.
Anyway, here is another way to do it, possibly in more transparent notation.
Begin by defining commutators and anti-commutators of $p$ and $q$ as simple superoperators,
$$
Q_+({\hat u}) = \frac{1}{2}\lbrace q {\hat u} + {\hat u} q \rbrace\\
Q_-({\hat u}) = q {\hat u} - {\hat u} q\\
\Pi_+({\hat u}) = \frac{1}{2}\lbrace p {\hat u} + {\hat u} p \rbrace\\
\Pi_-({\hat u}) = p {\hat u} - {\hat u} p
$$
where ${\hat u}$ is some arbitrary operator. Then observe that
1) $Q_\pm$ and $P_\pm$ satisfy very nice commutation relations of their own:
$$
\left[ Q_-, \; Q_+\right] = \left[ \Pi_-, \; \Pi_+\right] = \left[ Q_-, \; \Pi_+\right] = \left[ Q_+, \; \Pi_- \right] = 0\\
\left[ Q_\pm, \; \Pi_\pm\right] = i
$$
2) In terms of $Q_\pm$ and $P_\pm$, your generator simplifies to
$$
{\mathcal L} = {\mathcal L}_0 + {\mathcal L}_1(t)\\
{\mathcal L}_0 = \Pi_+ \Pi_- + \omega^2 Q_+ Q_- \\
{\mathcal L}_1(t) = -A(t) Q_-^2 + C(t) \Pi_+ Q_-
$$
Note: There is nothing new here, and in principle one may dispense entirely with points (1) and (2), but the compact notation and the simple algebra do make superoperator calculus much more straightforward.
3) As anticipated, it is now easy to see that ${\mathcal L}_0$ and ${\mathcal L}_1$ also have simple nested commutation relations:
$$
\left[ {\mathcal L}_0, \;{\mathcal L}_1(t)\right] = 2iA(t) \Pi_+ Q_- - iC(t)\left(\Pi_+^2 - \omega Q_-^2 \right) \equiv {\mathcal K}_1(t)\\
\left[ {\mathcal L}_0, \left[ {\mathcal L}_0, \;{\mathcal L}_1(t) \right] \right] = \left[ {\mathcal L}_0, \;{\mathcal K}_1(t)\right] = 4 \left[ A(t) \left(\Pi_+^2 - \omega Q_-^2 \right) + \omega^2 C(t)\; \Pi_+ Q_- \right] \equiv {\mathcal K}_2(t)\\
\left[ {\mathcal L}_0, \left[ {\mathcal L}_0, \left[ {\mathcal L}_0, \;{\mathcal L}_1(t) \right] \right] \right] = \left[ {\mathcal L}_0, \;{\mathcal K}_2(t)\right] = 2\omega^2 {\mathcal K}_1(t)\;,\;\;\;\text{etc.}
$$
where
$$
\left[ {\mathcal K}_1(t), \; {\mathcal K}_2(t)\right] = 0\;,
$$
and also
$$
\left[ {\mathcal L}_1(t), \left[ {\mathcal L}_0, \;{\mathcal L}_1(t) \right] \right] = \left[ {\mathcal L}_1(t), \;{\mathcal K}_1(t)\right] = 0\\
\left[ {\mathcal L}_1(t), \left[ {\mathcal L}_0, \left[ {\mathcal L}_0, \;{\mathcal L}_1(t) \right] \right] \right] = \left[ {\mathcal L}_1(t), \;{\mathcal K}_2(t)\right] = 0\;, \;\;\;\text{etc.}
$$
If $A$ and $C$ were time-independent, one could apply the Zassenhaus formula to the evolution $e^{{\mathcal L} t} = e^{\left({\mathcal L}_0 + {\mathcal L}_1\right) t}$,
$$
e^{\left({\mathcal L}_0 + {\mathcal L}_1\right) t} = e^{{\mathcal L}_0 t} e^{{\mathcal L}_1 t} e^{-(t^2/2!) \left[ {\mathcal L}_0, \;{\mathcal L}_1\right]} e^{(t^3/3!) \left(\left[ {\mathcal L}_0, \left[ {\mathcal L}_0, \;{\mathcal L}_1 \right] \right] + 2 \left[ {\mathcal L}_1, \left[ {\mathcal L}_0, \;{\mathcal L}_1 \right] \right] \right)} \cdot \dots
$$
and observe that all factors after the first two, $ e^{{\mathcal L}_0 t} e^{{\mathcal L}_1 t}$, mutually commute and that the exponents reduce to a single term each and can be eventually summed up.
On the other hand, the time-dependence of $A$ and $C$ turns the overall evolution into
$$
e^{{\mathcal L} t} = {\mathcal T} e^{\int_0^t{d\tau\;\left({\mathcal L}_0 + {\mathcal L}_1 \right)} }
$$
where $\mathcal T$ stands for time ordering, and the procedure complicates as usual. But perhaps in simpler cases it may still provide some good guidelines. For the current problem the final target is of course an explicit algebraic expression for the time-dependence of the number operator.
