How to model Phase noise?
Does this example prove to you that the spectrum spreads ?
The spectral plot uses dB (I find that more convenient) on the Y-axis, so consider it a logarithmic scale.

You can see that I need to add A LOT of noise to the phase.
I think this has to do with the fact that the noise I add is white noise.
In practical solutions this noise usually has a 1/f character near the carrier.

I used Octave for this as I do not have Matlab available (and Octave works for me and is free :-) )

clear all;clc;

Fin =1;
Ttotal = 1000;
nop = 100000;

Tsample = Ttotal/nop;
t = 0:Tsample:Ttotal;
Fsample = nop/Ttotal;

x  = sin(2*pi*Fin*t); 
xn  = sin(2*pi*Fin*(t + 0.1*randn(size(t))));

subplot(2, 1, 1);
points = 500;
plot(Fsample*t(1:points),xn(1:points),'r',Fsample*t(1:points),x(1:points),'b')
NFFT = 2^nextpow2(nop); % Next power of 2 from length of y
SX  = fft(x,NFFT)/nop;
SXn = fft(xn,NFFT)/nop;
dBSX = 20 * log( abs(SX));
dBSXn = 20 * log( abs(SXn));
f = Fsample/2*linspace(0,1,NFFT/2+1);
subplot(2, 1, 2);
plot(f,dBSXn(1:NFFT/2+1),'r',f,dBSX(1:NFFT/2+1),'b')
axis([0.5,1.5,-200,0]);