Gaussian function power spectrum (Matlab)
The problem is your calculation of the theoretical Fourier values.

When you do a DFT with N points, the interval between points on the frequency axis is  \$ \Delta f = 1/(N t_s) \$ where \$ t_s\$ denotes the sampling time. 

Therefore you have to calculate the theoretical Fourier values with
$$
\sqrt{\frac{\pi}{a}}\exp{-(\frac{(\pi k)^2}{(Nt_s)^2\cdot a})}
$$

Use nu=[-n:1:n]; and y=(sqrt(pi/a))*exp((-pi^2*nu.^2)/((2*n+1)^2*a)); in your code.