Before anything, let me just rewrite your TF from linear frequency ($f$, in Hertz) to angular frequency ($\Omega = 2\pi f$, in rad/sec). The notation $\Omega$ is usually adopted in the field of DSP to differ it from the discrete angular frequency ($\omega = \Omega T_s$, in rad), where $T_s$ (in seconds) is the sampling period. Therefore, the TF is given by $$ G(j\Omega) = e^{- \frac{(\sigma \Omega)^2}{2}} $$
To plot the DTFT and the DFT, you must first define the parameters of your problem. The Gaussian filter, as you stated, has $\sigma$ as parameter, which is the standard deviation. By setting this to, say, $\sigma = 0.1$, we can plot this impulsive response and this FT.
FT
This is a advisable step as it is sensible to take a look at this spectrum before choosing a sampling rate in order to avoid aliasing. In Matlab, you can express a Dynamic System Model in tons of manners. In your case, you already have the analytical expression $G(j \Omega)$, so lets make use of it:
%% Continuous-time plot % parameters Np1 = 401; % number of points for the plot of g(t) and |G(f)|,i.e., a total of N+1 samples t = linspace(-2,2,Np1); % time base (from -2s to 2s) Ts = t(2)-t(1); % sampling period (10 ms) fs = 1/Ts; % sampling rate (100 Hz) Omega = linspace(-pi/Ts,pi/Ts,Np1); % frequency base (from -100π rad/s to 100π rad/s) sigma=0.1; variance=sigma^2; g_t = 1/(sqrt(2*pi*variance))*(exp(-t.^2/(2*variance))); % g(t) G_Omega = exp(-(sigma*Omega).^2 / 2); % G(jΩ) % plot subplot(2,1,1) plot(t,g_t,'b'); grid on; title('g(t)'); xlabel('Time(s)'); ylabel('Amplitude'); subplot(2,1,2) plot(Omega, abs(G_Omega),'r'); grid on; title('|G(j\Omega)|'); xlabel('Frequency (rad/sec)') ylabel('|G(j\Omega)|'); xlim([Omega(1), Omega(end)]);

You might want to plot $|G(j\Omega)|$ in semilogarithm scale. In this case, assuming that $g(t)$ express a root-power quantity, such as a voltage or current in a electronic circuit, $|G(j \Omega)|$ in dB is given by: $$ |G_{dB}(j \Omega)| = 20 \log(|G(j \Omega)|) $$
Anyway, one can see that that most of the effective bandwidth of $G(j\Omega)$ is something between $-30$ and $30$ rad/sec. By denoting the one-sided bandwidth of $G(f)$ as $\Omega_B = 30$ rad/sec, from the Nyquist criterion, you must select a sampling rate such that $$ \Omega_s = \frac{2\pi}{T_s} \geq 2\Omega_B $$ By adopting $f_s = 1/T_s = 100 Hz \therefore \Omega_s = 2\pi 100 $ rad/s, we can fulfill the Nyquist requirement.
DTFT
One way to implement a digital filter is by sampling the impulse response, $g(t)$, in order the obtain the same impulse response of the equivalent analog filter. Such an approach is usually called impulse invariance method. However, we usually normalize the signal samples by $T_s$ so that we cancel the $T_s$ is the DTFT and get the very same transfer function of the analog filter (vide Oppenheim, eq (4.20)). Hence, by defining $$ g[n] = T_s g(nT_s) $$ its DTFT is given by $$ G(e^{j\omega}) = \sum_{k = - \infty}^{\infty} G\left[j\left(\frac{\omega}{T_s} - \frac{2\pi k}{T_s} \right)\right], $$
In the particular case of no aliasing, we have that
$$ \begin{align} G (e^{j\omega}) = G\left( j\frac{\omega}{T_s}\right) & \textrm{ , for }|\omega| \leq \pi \end{align} $$
It is obvious that $G_s (e^{j\omega})$ is actually defined in $(-\infty, \infty)$, but we usually represent this signal only in $[-\pi, \pi]$ as it is periodic in $2\pi$.
%% DTFT G_DTFT = G_Omega; omega = Ts .* Omega; plot(omega, abs(G_DTFT), 'b'); title('|G(e^{j\omega})|'); xlabel('\omega (in rad)'); ylabel('Amplitude'); grid on; xlim([-pi,pi]); xticks([-pi,-pi/2, 0, pi/2, pi]); xticklabels(["-\pi","-\pi/2", "0", "\pi/2", "\pi"]);

As you can see, plotting the DTFT is ridiculously easy (as long as you ensure that aliasing does not occur!) since it just involves frequency scaling (if you don't normalize the samples by $T_s$, it involves a scaling in the amplitude as well).
DFT
For the DFT, you actually have two way outs. For the first one, note that (vide Oppenhneim, again)
$$ \left.G[k] = G(e^{j\omega})\right\vert_{\omega = \frac{2\pi k}{N+1}} $$ where $N+1$ is the number of samples of the sequence $\{G[k]\}$, which is also equal to the number of samples of the sequence $\{g[n]\}$. Since what you have on computer is samples of $G(e^{j\omega})$ instead of itself, strictly speaking, you already made it. The plot() function makes the illusion, by interpolating the samples, that you have a continuous curve instead of a vector of values. Since we have plotted $G(e^{j\omega})$ from $-\pi$ to $\pi$, let us define $k \in \left\{-\frac{N}{2}, -\frac{N}{2}+1, \cdots, \frac{N}{2}-1, \frac{N}{2}\right\}$.
%% DFT 1th method - by hand G_DFT = G_DTFT; k = -(Np1-1)/2:(Np1-1)/2; stem(k, abs(G_DFT), 'b'); title('|G[k]|'); xlabel('k'); ylabel('Amplitude'); xticks([-(Np1-1)/2, -(Np1-1)/4, 0, (Np1-1)/4, (Np1-1)/2]); xticklabels(["-N/2", "-N/4", "0", "N/4" "N/2"]); grid on;

The second way out, which kinds of serves as a comparison with the Matlab built-in function, is directly applying the FFT on the samples $\{g[n]\}$ (remember that the FFT is a efficient implementation of the DFT computation, and the final result is exactly the same). The computation of the FFT is given by: $$ G[k] = \frac{1}{N+1} \sum_{n = 0}^{N} g[n] W_{N+1}^{nk}, $$ where $W_{N+1} = e^{-j\frac{2\pi}{N+1}}$ is called twiddle factor.
The magic is that, by computing $G[k]$ from the previous equation, it shall exactly be the samples of $X(e^{j\omega})$ in the points $\omega = \frac{2\pi k}{N}$ as the relationship between DTFT and DFT must match. Just pay attention that your vector of samples must begin with $g[0]$. Likewise, the sequence $\{G[k]\}$ returned by fft() also begins with $G[0]$. Shifting it in order to obtain the very same plot is easily done by fftshift().
%% DFT 2th method - by using the built-in FFT function g_k = Ts * g_t; % g[k] G_k = fft(fftshift(g_k)); % G[k] (the DFT of g[k]) stem(k, abs(fftshift(G_k)), 'b'); title('|G[k]|'); xlabel('k'); ylabel('Amplitude'); xticks([-(Np1-1)/2, -(Np1-1)/4, 0, (Np1-1)/4, (Np1-1)/2]); xticklabels(["-N/2", "-N/4", "0", "N/4" "N/2"]); ylim([0, 1]); grid on;
