Matlab: Correct method of generating Rayleigh fading channel
Question 1: 

Yes, it is a multipath Rayleigh Channel with $ 1000 $ taps. You need to know that when you generate complex coefficients, and the real part and the imaginary part of this coefficients are randomly chosen according to a Gaussian distribution, the magnitude will have a Rayleigh distribution and the phase will have a uniform distribution between $ 0 $ and $ 2 \pi $. The problem with this method, is that you don't have control on how the channel varies: the channel in a given instant of time will have no correlation with the channel in the previous instant. To add this correlation, we use the others methods as Clark's and Jake's models.

Question 2:

Look what he wrote in his website, in the assumptions section:


  
  The channel is flat fading – In simple terms, it means that the multipath channel has only one tap. So, the convolution operation reduces to a simple multiplication.
  


In his code, he is performing a multiplication, not a convolution, when applying the Rayleigh Channel.

He generates $ N $ coefficients, but what he is trying to do is simulate the channel variation in time, so each symbols is multiplied by an attenuation (inside the vector h). The vector h is not a vector of a filter coefficients, but a vector that stores the channel attenuation for each symbol.