Scaling problems with IFFT in Matlab

I'm studying the IFFT in Matlab by applying it to a Gaussian. According to Wikipedia tables, the Fourier transform pair would be

F(w) = sqrt(pi/a) * exp(-w^2/(4a))

in frequency, and

f(t) = exp(-at^2)

in time. I modified the code in a previous question plus Cris Luengo's answer to perform this IFFT.

a = 0.333;

ts = 1e4; % time sampling
L = 1000*ts; % no. sample points
ds = 1/ts;
f = -floor(L/2):floor((L-1)/2); % freq vector
f = f/ts;

w = 2*pi*f; % angular freq
Y = sqrt(pi/a)*exp(-w.^2/(4*a));

y = ts*ifftshift(ifft(fftshift(Y)));

t = (-L/2:L/2-1)*ts/L; % time vector

f = exp(-a*t.^2); % analytical solution

figure; subplot(1,2,1); hold on
plot(t,real(y),'.--')
plot(t,real(f),'-')
xlabel('time, t')
title('real')
legend('numerical','analytic')
xlim([-5,5])
subplot(1,2,2); hold on
plot(w,imag(y),'.--')
plot(w,imag(f),'-')
xlabel('time, t')
title('imag')
legend('numerical','analytic')
xlim([-5,5])

When I compare the result of IFFT with the analytical expression, they don't seem to agree:

I'm not sure where the mistake is. Have I scaled the IFFT properly? Is there an error in how I define the linear/angular frequency?

Edit: For some reason, when I define L=ts^2 the analytical and numerical solutions seem to agree (L = no. sampling points, ts = time sample).

Starting from the analytical solution, let's rephrase things a bit. You have a sampling of the function f(t) = exp(-a*t^2), and the way you've constructed the analytical answer you're collecting L=1000*ts=1e7 samples at a sampling rate of Ts=ts/L=1e-3. This means that your sampling frequency is Fs=1/Ts=1e3.

Since you want to compare against results obtained with fft/ifft, you should be considering digital or discrete frequencies, meaning the values you define for your transform will correspond to the digital frequencies

frd = (-L/2:L/2-1)/L;

Mapping this to angular frequencies, we have:

w = 2*pi*frd;

But when you're trying to compute the values, you also need to keep in mind that these frequencies should represent samples of the continuous time spectrum you're expecting. So you scale these values by your sampling frequency:

Y = sqrt(pi/a)*exp(-(Fs*w).^2/(4*a));
y = Fs*ifftshift(ifft(fftshift(Y)));

When you compare the analytical and computed answers, they now match.

The short answer to your question, given this, is that you are scaling y incorrectly at the end. You're scaling it by ts, which is 1e4, but you need to be scaling it by the sampling frequency which is Fs=1e3. That's why you end up off by a factor of 10.