Tacking on a term $e^{-iEt/\hbar}$ is the correct interpretation of the Schrödinger equation $$i\hbar |\partial_t \Psi\rangle = \hat H |\Psi\rangle$$only for those eigenstates for which $$\hat H |\Psi\rangle = E|\Psi\rangle,$$as otherwise you do not know what value of $E$ should be used to substitute. Hypothetically you can still do it, but you pay a very painful cost that the $E$ is in fact a full-fledged operator and you therefore need to exponentiate an operator, which is nontrivial.
If this is all sounding a bit complicated, please remember that QM is just linear algebra in funny hats, and so you could get an intuition for similar systems by just using some matrices and vectors, for example looking at $$i\hbar \begin{bmatrix} f'(t) \\ g'(t) \end{bmatrix} = \epsilon \begin{bmatrix} 0&1\\1&0\end{bmatrix} \begin{bmatrix} f(t) \\ g(t)\end{bmatrix}.$$One can in fact express this as $$\begin{bmatrix}f(t)\\g(t)\end{bmatrix} = e^{-i\hat H t/\hbar} \begin{bmatrix} f_0\\ g_0\end{bmatrix},$$ but one has to exponentiate this matrix. That is not hard because it squares to the identity matrix, causing a simple expansion, $$\begin{bmatrix}f(t)\\g(t)\end{bmatrix} = \cos(\epsilon t/\hbar) \begin{bmatrix} f_0\\ g_0\end{bmatrix} - i \sin(\epsilon t/\hbar) \begin{bmatrix} g_0\\ f_0\end{bmatrix}. $$ One can then confirm that indeed this satisfies the Schrödinger equation given above. One can also immediately see that this does not directly have the form $e^{-i\epsilon t/\hbar} [f_0; g_0],$ but how could it? That would be a different Hamiltonian $\hat H = \epsilon I.$
But, with some creativity, one can see that if $f_0 = g_0$ those two remaining vectors would be parallel, or if $f_0 = -g_0$, and one can indeed rewrite this solution in terms of those eigenvectors of the original $\hat H$ as $$\begin{bmatrix}f(t)\\g(t)\end{bmatrix} = e^{-i\epsilon t/\hbar} \alpha \begin{bmatrix} 1\\ 1\end{bmatrix} + e^{i\epsilon t/\hbar} \beta \begin{bmatrix} -1\\ 1\end{bmatrix}. $$ So the trick to more easily finding general solutions is to find these eigenvectors first and then form a general linear combination of those eigenvectors once they have been multiplied individually by their time dependence. Then for a given initial state, we need to find the $\alpha$ and $\beta$ terms: in this case it is simple enough by looking at $t=0$ where $\alpha - \beta = f_0$ while $\alpha + \beta = g_0.$
Similarly for your Hamiltonian $\hat H = \hat p^2/(2m) = -\frac{\hbar^2}{2m}\frac{\partial^2~}{\partial x^2},$ you know that the eigenvectors are plane waves, $$\phi_k(x) = e^{ikx}.$$You know that you can then add time dependence to them in the obvious way, $$\Phi_k(x, t) = e^{i(k x - \omega_k t)},$$ where of course $$\hbar \omega_k = \frac{\hbar^2k^2}{2m}.$$ So the eigenvector story is just beautifully simple for you to do, all you need is the ability to take derivatives of exponentials.
The part of the story that is more complicated is assembling an arbitrary $\psi(x)$ as a sum of these exponentials. However while it is complicated it is not impossible: you know from Fourier's theorem that $$\psi(x) = \frac{1}{2\pi}\int_{-\infty}^{\infty} dk ~e^{i k x} \int_{-\infty}^\infty d\xi ~e^{-ik\xi} ~\psi(\xi).$$ Let your eyes glaze over the second integral and see it as just what it is, some $\psi[k]$ coefficent in $k$-space. What we have here then is a sum—a continuous sum, but still a sum!—of coefficients times eigenfunctions:$$\psi(x) = \int_{-\infty}^{\infty}\frac{dk~\psi[k]}{2\pi}~\phi_k(x).$$
And we know how to Schrödinger-ize such a sum, we just add $e^{-i\omega_k t}$ terms to each of the eigenfunctions, turning $\phi_k$ into $\Phi_k.$ So we get,
$$\Psi(x, t) = \frac{1}{2\pi}\int_{-\infty}^{\infty} dk ~e^{i (k x - \omega_k t)} ~\psi[k].$$
You do not have to do it this way, you can try to do some sort of $$\exp\left[-i \frac{\hbar t}{2m} \frac{\partial^2~}{\partial x^2}\right] e^{-a x^2}$$
monstrosity, expanding the operator in a power series and then seeing whether there are patterns you can use among the $n^\text{th}$ derivatives of Gaussians to simplify. But the operator expansion way looks really pretty difficult, while the eigenvector way is really easy.
The reason it is really easy is that both $\hat H$ and $i\hbar \partial_t$ are linear operators: they distribute over sums. So if you are still feeling queasy about this procedure, convince yourself by just writing it out: calculate this value $$0 = \left(i\hbar \frac{\partial~}{\partial t} + \frac{\hbar^2}{2m}\frac{\partial^2~}{\partial x^2}\right) \frac{1}{2\pi} \int_{-\infty}^\infty dk~\psi[k] ~e^{i (k x - \omega_k t)}.$$ Notice that it holds with pretty much no restriction on the actual form of $\psi[k]$ so that you only need to choose coefficients $\psi[k]$ such that $\Psi(x, 0) = \psi(x).$
