AWGN and Rayleigh fading in slow flat fading channels in MATLAB
As I described in this answer, it is generally assumed in the Rayleigh fading model that the energy is conserved on average (i.e. that the fading coefficient $R$ is such that $E\left\{R^2\right\} = 1$). As a result you would need to generate the Rayleigh distribution from two normal distribution with standard deviation of $1/\sqrt{2}$ (instead of with unit standard deviation). This can be achieved with:

rayleigh = sqrt( (random('norm',0,1/sqrt(2)))^2 + (random('norm',0,1/sqrt(2)))^2 );


or alternatively, since you are already using random from the statistics toolbox, with:

rayleigh = random('Rayleigh',1/sqrt(2));


Other than that, the expressions you have properly represent a slow flat Rayleigh fading channel for the transmission of a one-dimensional signal (such as BPSK, M-PAM, ...)