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).$$$$\Psi(q,t) = e^{-iHt/\hbar}\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.