This will get you part of the way. You need to figure out the reason your code does not reproduce the published plots. This aspect of communications engineering is not in an area of my expertise, so I cannot help you with those details.
The (Slightly Changed) Code:
fc = 553e6;
w = 0:2*pi*fc*1e-2:100*2*pi*fc;
pau = 0.33;
tau = [0.2 0.5 1.0 2.0]*1e-9; % in nano second
[W,Tau] = meshgrid(w, tau);
% beta_square = zeros(1,length(w));
beta_square = @(w,tau) 1 + 2*pau*cos(w.*tau) + pau^2;
theta = @(w,tau) atand(-pau*sin(w.*tau)./(1+pau*cos(w.*tau)));
figure,plot(w,10*log10(beta_square(W,Tau)))
set(gca, 'XLim',[4 6]*1E+10)
figure,plot(w,theta(W,Tau))
set(gca, 'XLim',[4 6]*1E+10)
Note that ‘beta_square’ is in units of decibels, so I added that (assuming that it represents a power, if it represents an amplitude instead, 10*log10() becomes 20*log10()). The phase is in degrees, so I also changed atan to atand.
The code works. You may have to review the paper the equations and published figures come from to determine the reason your code is not reproducing those figures.
