One alternative way to go about it would be to use the already known static solutions of the harmonic oscillator $\phi_n(q)$. They are eigenfunctions, so they evolve in time by picking up a phase factor $e^{i E_n t/\hbar}$, with $E_n = \hbar \Omega (n + 1/2)$. These functions are a complete basis of the Hilbert space, so you can especially decompose the initial conditions as
$$ \Psi_0(q) = \Psi(q,t=0) = \sum_{n\ge 0} \langle\phi_n |\Psi_0\rangle \phi_n(q).$$
You can thus obtain the solution as
$$\Psi(q,t) = e^{iHt}\Psi_0(q) = \sum_{n\ge 0} e^{iE_nt/\hbar} \langle\phi_n |\Psi_0\rangle \phi_n(q).$$
This shifts the task from getting a new solution of the equation to just calculating $\langle \phi_n | \Psi_0\rangle$, which is feasible with Gaussian integrals and using a series expansion for the Hermite polynomials (which should converge).

There isn't much elegance in this solution. It's quite technical and maybe the detour over the Riccati equation yields a more useful result.

The ansatz you use leads to a Riccatti equation, but the anstaz has a flaw: it is not normalized. The norm is time-dependent, 
$$ \langle \Psi(t) | \Psi(t)\rangle = |A|^2 \sqrt{\frac{\pi}{\lambda(t)}}.$$
This can not yield the correct solution because the Schrödinger equation is unitary, that is the norm of $\Psi$ is necessarily time-independent.