From the top of my head:
Fc=1e6; % carrier 1MHz
Fm=1e3; % single tone - 1kHz
fs=2.2*Fc; % sampling frequency (Nyquist limitation!)
t=0:1/fs:1;
Y=zeros(1,length(t));
for n=-20:20
Y=Y+besselj(n, B).*cos( 2*pi*(Fc-n*Fm).*t);
end
Should work (right now I can't vertify this code - lack of Matlab)