Rephrasing OP's Post for Clarity
It wasn't stated, but I assume the OP is attempting to create a one-tap channel filter with:
h=(1/sqrt(2))*(randn(N,1)+1i*randn(N,1))
as an array of $N$ samples,
and a four-tap channel filter with (I renamed for clarity)
g=(1/sqrt(2))*(randn(N,4)+1i*randn(N,4))
as an array of N 4 sample vectors.
The OP uses $N=1000$, as each sample in $h[n]$ or 4 samples in vector $g[n]$, where each sample is the the sum of two Gaussian random variables in quadrature representing a Rayleigh model: it's magnitude is Rayleigh distributed and its phase is uniformly distributed. With the one-tap solution, each $h[n]$ is a single complex value, and with the four-tap solution, each $g[n]$ are four complex values, as the coefficients of an FIR filter.
It is not clear what the OP really intends to do with the $1000$ samples - this ultimately represents 1000 different experiments in transmitting from the transmitter to the receiver, and how often we change $h[n]$ or $g[n]$ depends on if we are simulating slow fading or fast fading environment, or if we just want to see the statistical distribution of all possible fades in either case (essentially running an experiment $N$ times, the system is ergodic). Which is used would be a factor if using this model to test something like equalizer convergence and how other temporal operations in our receiver might behave in the presence of such time varying fading (such as clock and carrier recovery, and automatic gain control).
Not covered in this post is Doppler spread which is also typically included in channel simulations for mobile applications or if objects as sources of reflection are moving relative to the transmitter and receiver.
Question 1
Does this create a frequency selective channel solely because I have
more taps?
Perhaps but possible 4 taps is insufficient. A flat fading channel is a single-tap FIR filter: we simply modify the gain and phase. If it is a "Rayleigh Channel" then we select from a zero-mean complex Gaussian random process a magnitude and phase and multiply the waveform by that. (How often we change from one $h[n]$ to the next determines fast or slow fading which is another parameter).
A frequency selective channel occurs when the delay spread of the channel (the time dispersion between the multiple paths which itself ends up giving us the randomly selected magnitude and phase we are using) exceeds the symbol duration. We can use an FIR filter to represent this spread occurring over multiple symbols given the sampling rate of the waveform. If the OP's waveform was sampled at four samples per symbol, the achieved delay spread would not exceed a symbol duration and therefore would not be frequency selective fading. The delay spread to use is based on the known statistics for the environment (see the background section below for further details.).
Question 2
In non-of these definitions we are taking into account the delay and
magnitude of each path
The delays are the coefficient position of the 4 tap filter. If we wanted to simulate $T=[2,3,4,5]$ then this would be a 6 tap filter with coefficients assigned from $h[n]$ (in the 1000,4 shape case) as $[0,0,h[0],h[1],h[2],h[3]]$. The first two zeros do nothing to effect the model, (added time delay with no additional distortion) so really can be omitted. The coefficients change both the magnitude and phase, as we desire.
Details on how to use OP's Samples in a Channel Simulation
It really didn't matter to generate $Nx1$ or $Nx4$ samples, all the samples given are independent and identically distributed samples from a complex Gaussian process. How many samples we need depends on what we are actually simulating in terms of frequency selective fading (how many samples to populate the filter) and fast vs slow fading (how often we change the samples in the filter. I will however continue to use $h[n]$ to refer to individual complex samples that may change with time for slow or fast fading and $g[n]$ to refer to a vector of samples representing the state of the channel filter that also may change with time depending on the rate of fading. The length of $g[n]$ however depends on the degree of frequency selective fading.
To account for frequency selective fading, simply populate an FIR filter where the time span of the filter in seconds (number of filter taps divided by sample rate) is equal to the delay spread in seconds desired (this would depend on environment such as urban, indoor, etc and frequency band). If there is no frequency selective fading (delay spread < symbol period) this will result in a single-tap filter.
How often we change the samples accounts for modelling fast or slow fading. Slow fading means the the samples are changed much slower than our symbol rate, and fast fading results when the samples are changed closer to the symbol rate. One solution for simulating fast vs slow fading is to perform a moving average on the generated fading samples $h[n]$ or each value in $g[n]$ to create a new samples with some coherence from sample to sample consistent with the rate of change of fading.
This approach will also give a "knob" to change between fast and slow fading by changing the number of samples $M$ in the moving average.
So in summary there is an FIR filter used that the waveform goes through to simulate a channel, this is a single tap filter for flat fading or a multi-tap filter for frequency selective fading. Each coefficient is selected as a sample of a complex Gaussian random process. The coefficients can be generated from a moving average to create a time varying filter to simulate fast or slow fading.
Do not be confused by the use of the FIR filter as the sum of delayed paths itself and the terminology of Rayleigh fading simulating multiple paths: Each tap of the filter already represents multiple paths of a Rayleigh channel within that coherence time (so a single tap solution is still many reflection paths, they just are not spread significantly enough in time to require a multiple tap filter). The fact that we use a Gaussian random variable to set the coefficient of the tap has already accounted for the summation of multiple random paths.
More Background
The following provides more details for those interested in further understanding Rayleigh Fading with regards to Frequency Selective or Flat Fading:
Rayleigh Fading
The distribution of the amplitude of a complex Gaussian zero mean signal is a Rayleigh distribution as you describe, this is because a complex Gaussian signal is the sum of two independent equally distributed Gaussian signals given as $I + jQ$ where I is the in-phase (Real) and Q is the quadrature (Complex).
The reason the received signal strength for certain fading channels are modeled as a Rayleigh distribution is due to the central limit theorem: As you sum independent random variables the resulting distribution approaches a Gaussian distribution.
In the case of fading channels that have received signal strength that can be modeled as Rayleigh, these would be channels with multiple independent paths with no direct path that would be significantly stronger from all the reflected paths. If the delay in each path is long enough (delay spread), each path would be uncorrelated from the other paths as received, and when combined (if enough paths) would appear as a complex zero-mean Gaussian random variable. The amplitude of this (our received strength) is Rayleigh distributed as described earlier.
In contrast is the fading environment where there is a dominant stronger direct path. The result of this is a complex Gaussian random variable with a mean (the direct path), and the amplitude of such a random variable follows a Ricean distribution.
Frequency Selective or Flat Fading
To be frequency selective fading or flat fading is independent of the distribution of the signal strength being Rayleigh. That only has to do with delay spread. (See figure below).
The Delay Spread of the channel is what determines the frequency selectivity of the fading (whether it is Rayleigh or not). Typical delay spreads for the macrocellular environment are shown in the figure below:
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$