Your approach will not model the actual channel but can be used to create the received waveform as having passed through a Rayleigh fading channel. I describe considerations to actually create the received waveform using your approach in either flat and frequency selective conditions below, which should give insight into why it would not be a channel simulation.
To simulate a channel
The Rayleigh channel can be simulated with an FIR filter with a significant number of taps that match the sample rate of the waveform and enough taps that the span of the filter exceeds the delay spread of the channel. As for "correct ways" of simulating a Rayleigh fading channel itself as a FIR filter, please refer to "Jakes' Channel Model" which uses a "sum of sinusoids" approach to modeling the Rayleigh channel. More information on Jakes' model and other similar approaches is available here.
Details on how to proceed with the "Dictionary Approach"
The dictionary would represent the effect of a "flat fading" Rayleigh channel as a "1 tap filter" as you described if your dictionary was generated at your symbol rate and you multiply (not filter!) each transmitted symbol with a sample from your dictionary. In this fashion each symbol would have an independent amplitude and phase from the Rayleigh channel and that fade would be constant over a symbol duration consistent with flat fading. However the fade would be completely independent from one symbol to the next, so in this case would poorly represent an actual channel.
One solution to that problem is to perform a moving average on the generated dictionary to create a new dictionary with some coherence from symbol to symbol. Much better would be to create an upsampled dictionary by $M$ and then do a moving average over $M$ samples. The higher $M$ is, the more accurately you will represent a true Rayleigh channel. In this approach the transmitted symbols would also need to be interpolated to match the rate in the dictionary. Your transmitted symbols would be interpolated for pulse shaping and other purposes in your modelling, so the final dictionary using this last approach described can be decimated at the output of the moving average filter to match the number of samples required in your waveform for all other reasons to accurately represent your model.
This approach will also give you the "knob" you are looking for to change between frequency selective and flat fading; what is described above by creating a symbol length moving average over $M$ samples, with $M$ samples per symbol, is at the threshold position between flat fading and frequency selective fading. If you decrease the span of the moving average to be less than the symbol rate, it will become frequency selective fading, as you will have uncorrelated fading samples within your symbol period which is frequency selective fading.
So in summary to create a received waveform that has passed through a Raylegih channel do the following:
1) Create your transmitted waveform (from $N$ complex baseband symbols) that is sufficiently upsampled to model your channel bandwidth of interest (usually for pulse shaping and spectral regrowth or whatever else you may be modelling; this will typically be between 2 and 10 samples per symbol). I will use $K$ to be the total number of samples, and lower case $m$ to be the upsample factor, so then the total number of samples will be $K = Nm$.
2) Create a vector of $KM$ samples of independent complex Gaussian values using the process you described, where $M$ is an additional smoothing factor.
3) Do a moving average over $L$ samples on the vector created in step 2 above. Choose $L$ based on modeling either flat fading or frequency selective fading where the threshold boundary between the two cases is when the moving average span is equal to the symbol duration ($L \ge mM$ is flat fading and $L < mM$ is frequency selective fading).
4) Decimate the $KM$ samples a the output of the moving average smoothing filter to match the sample rate of your $K$ interpolated transmitted waveform to create the final dictionary. (Decimate by $M$).
5) Finally multiply each of the $K$ samples in the transmitted waveform by each of the $K$ samples in the final dictionary as a sample by sample multiplication to create the received waveform as having gone through a Rayleigh fading channel.
The higher $M$ is, the more accurate the channel result would be but I would think a factor between 6 to 10 would both be easy to manage and be sufficient for creating the coherence from sample to sample. For more frequency selectivity, an even higher $M$ may be needed using this approach. Ultimately with high frequency selectivity each sample is decorrelated and the use of the moving average is unnecessary (the degree of frequency selectivity would be controlled by the interpolation factor of your waveform). I do like the use of the moving average as it gives you the "knob" to control the degree of fading between flat and frequency selective.
To be most accurate (although a lot more complicated), you could optionally vary L over the course of the moving average process within the range consistent that is consistent with your desired fading environment (flat or frequency selective)...this changes the coherence itself over time which would more accurately represent channel conditions, but that may be a second order effect; meaning not critical to your model).
N=1; h=(1/sqrt(2))(randn(N,1)+1i*randn(N,1));. Could you explain what are $N$ and $M$ in your code? $\endgroup$