How can I shape my time channel response, in order to have Gaussian shaped Doppler Spread?
See  section 2.4.3 of this reference https://web.stanford.edu/~dntse/Chapters_PDF/Fundamentals_Wireless_Communication_chapter2.pdf

If the doppler spectrum has to be gaussian, the auto-correlation of tap gains should be gaussian (which is correctly mentioned in other answer but I somehow felt more details were not captured). For a coherence time $T_c$, proportional to doppler spread $1/D_s$, will be expressed as smallest amount of $n$ for which your tap at $l$ will be significantly different from tap at $l+n$. So for example, if $n=5$, gain of taps from $0$ till $4$ will be having a guassian distribution, while the gain at tap $0$ will be uncorrelated with tap at $5$.

If you model your taps as Rayleigh fading channel, your PSD will not be Gaussian as given in eq 2.60 of the reference. You may need to generate your tap gains such that their auto-correlation function follows gaussian distribution.

Since you are using MATLAB

raychan = comm.RayleighChannel('DopplerSpectrum',{doppler('Gaussian')});


This makes sure that each of the channel taps gain auto-correlation (and PSD) follows Gaussian distribution.