Hello,
I don't know mathlab, but maybe this may help you.
The Rayleigh distribution is the sum of two orthogonal Gaussian variables.
When you start with two carriers, 90 degrees out of phase. Modulate both carriers with uncorrelated Gaussian noise sources and add the outputs from the multipliers. This will give you a Rayleigh faded carrier.
You need to put low pass filtering at both outputs of the Gaussian noise sources to model the speed of the envelope and phase variations. A slow varying path requires a low corner frequency for the LPF, while a higher setting for the low pass filter enables fast varying fading (driving fast in a city).
When you add a constant to (one of) the Gaussian noise sources, before multiplying with the carriers, you create less deep fades (in fact there is a dominant propagation path in that case, I think a Rician fading case).
When you have a modulated carrier already, you have to generate the 90 degrees out of phase signal with a 90 degrees phase shifter. This gives you two modulated carriers that you can use as input for the multipliers.