You need to vectorise (do element-wise calculations) in your ‘modz’ calculation. This means substituting ‘./’ for /:
modz=sqrt((w.^0.5 + 20*cos(pi/8)*w.^0.25 + 100)./(1+ 2*cos(pi/8)*w.^0.25 +w.^0.5));
↑
That will work.
