It is hard to know what exact the professor has done, but from what I have understood made an effort to help the situation.
Q1. A Gaussian function is chosen to represent a freely expanding wave function for several reasons:
(i) The Gaussian function represents a normal probability distribution function. Since $|\psi(x)|^2$ represents a probability distribution function for a huge class of particles moving in a similar way, and having momentum within a certain range, the large numbers theorem points towards a Gaussian function, to represent the wave function of particles within a region of space determined by the width, $\sigma$, of the Gaussian function. The $\sigma$ also happens to be the standard deviation, i.e. the uncertainty in the position of the particle.
(ii) The expectation value of the position of the particle turns out to be the value of $b$ in your Gaussian wave packet.
(iii) The Gaussian wave packet also contains the wavy bit shown by the phase factor
$e^{ikx}$.
This is what makes Gaussian wave packets an excellent representation of particles, which are known to be located within some region of width $w$, but nevertheless they are all moving as plane waves while the packet travels along in space.
(iv) The wave packet is made of an infinitely large number of momentum values in a momentum width which relates to $\sigma$ by a Fourier transformation.
Q2. So to put some order in all these, let us consider the general Gaussian function
$\psi(x)=\psi_0e^{-A(x-x_0)^2+Bx+C}$
which has a Fourier transformation
$\psi(k)=\psi_0\sqrt{\frac{\pi}{A}}e^{\frac{B^2}{4A}+Bx_0+C}$.
The normalisation constant $\psi_0$ is given by an integration in the range [-$\infty, +\infty$] and has the value $\sqrt{\frac{A}{\pi}}$ but leave it as $\psi_0$.
Q3. To make contact with your professor, try to apply these results using the following:
$A=\frac{1}{2\sigma^2}$
$B=ik$
$x_0=b$
$C=0$
For question v1:
Since you're integrating the non-negative function $(x^2-l^2)^4$, you shouldn't get zero. Your mistake must be your expression for the antiderivative. Expanding out the integrand is probably a safe way to start.
There's something funky about the supposed answer for $N$. Look at the units, along with the entire expression for $\psi$. The units of $\psi$ should be such that $\left | \psi \right |^2 dx$ is dimensionless. That tells you what the dimension of $N$ ought to be given that $\psi\sim l^4$. The correct solution you gave is off.
Best Answer
Given a wavefunction, $\psi(x,0)$, the corresponding quantity, $\psi(x,t)$, is obtained from the integral: $$ \boxed{\psi(x,t) = \int_{-\infty}^{\infty}dy\,Z[x,t;y,0]\psi(y,0)\,} $$ where $Z[x,t;y,0]$ is a path integral, $\int_{x,0}^{y,t} \mathcal{D}x\,e^{\frac{i}{\hbar}I[x]}$, and $I[x]$ the single particle action of interest. For free non-relativistic particles, $I[x]=\int_0^td\tau\frac{m}{2}(\frac{dx}{d\tau})^2$, and the correctly normalised path integral reads: $$ Z[x,t;y,0] = \sqrt{\frac{m}{2\pi i\hbar t}}\,\exp\Big(i\frac{m(x-y)^2}{2\hbar t}\Big). $$ The initial wave packet is of the form: $$ \psi(y,0)=Ne^{\frac{i}{\hbar}ky}e^{-\frac{1}{2}b(y-x_0)^2}, $$ where you can read off the $b,N$ of interest from your formula, e.g. $b=1/(2\sigma_x^2)$.
The objective is to compute $\psi(x,t)$, which amounts to carrying out the $y$ integral in the above boxed formula. This is a Gaussian and so can easily be carried out. You'll primarily need the result: $$ \int_{-\infty}^{\infty} da\,e^{iza^2} = \sqrt{\frac{\pi i}{z}},\quad {\rm if }\quad {\rm Im}\,z\geq0,\quad z\neq0, $$ which follows from Cauchy's theorem applied to $\oint_C da\,e^{iza^2}=0$, where the contour $C$ is that in the figure:
It is convenient to write $z=\rho e^{i\phi}$ (with $0\leq \phi<\pi$), and then given any allowed $\phi$ you choose $\theta$ such that $\phi+2\theta=\pi/2$. This enables you to take $R\rightarrow \infty$ while ensuring the angular contour contributions vanish, and you are left with $\int_{-\infty}^{\infty}da\,e^{iza^2}=\sqrt{i}e^{-i\phi/2}\int_{-\infty}^{\infty}da\,e^{-\rho a^2}$. Evaluating the remaining Gaussian integral and rewriting the result in terms of $z=\rho e^{i\theta}$ leads to the displayed result.
The next step is to redefine your integration variable, $a\rightarrow a'=a-b$ (with $b\in \mathbb{R}$), $$ \int_{-\infty}^{\infty} da\,e^{iz(a+b)^2} = \sqrt{\frac{\pi i}{z}}, $$ which implies: $$ \int_{-\infty}^{\infty} da\,e^{iza^2+2izba} = \sqrt{\frac{\pi i}{z}}\,e^{-izb^2},\quad {\rm if }\quad {\rm Im}\,z\geq0,\quad z\neq0. $$ We are almost there. The only thing remaining to do is to notice that both the left- and right-hand sides of this last expression are analytic in $b$, so we can extend $b$ into the complex plane, $b\rightarrow w=b+ib'$. Finally, we redefine $w\rightarrow zw$. Therefore, for generic complex numbers $z$ (with ${\rm Im}\,z\geq0$, $z\neq0$) and $w$, $$ \boxed{\,\int_{-\infty}^{\infty} da\,e^{iza^2+2iwa} = \sqrt{\frac{\pi i}{z}}\,e^{-iw^2/z},\quad {\rm if }\quad z,w\in \mathbb{C},\quad {\rm and}\quad{\rm Im}\,z\geq0,\quad z\neq0\,\,\,} $$ This is all you need in order to evaluate the above integral that gives you $\psi(x,t)$. You just read off $z,w$ and $a$ as appropriate for your case of interest and apply this result. (Of course the same formula applies for general $x_0$.) The result is: $$ \boxed{\psi(x,t) = \frac{1}{\sqrt{1+i\frac{\hbar b }{m}t}}\exp\Bigg\{\!\!-\frac{i}{\hbar}\frac{k^2t}{2m}\frac{\big[1+i\frac{\hbar b}{k}(x-x_0)\big]^2}{1+i\frac{\hbar b}{m}\,t}\Bigg\}\psi(x,0)\,\,} $$ Recall that $ \psi(x,0)=Ne^{\frac{i}{\hbar}kx}e^{-\frac{1}{2}b(x-x_0)^2}. $