帮看程序哪里有问题，谢谢！

jifang29

clc
clear
close all
W=1.0e+004 *[ 0.1112 0.1492 0.1945 0.2472 0.3098 0.3830 0.4690 0.5947 0.6890 0.8294 0.9940 1.1878 1.4159 1.6836 1.9990 2.3700];%截止频率
b(1)=tan(pi/4-mod(W(1)/2,2*pi));
for i=1:length(W)
    if (i==1)
       y(i)=phase1(W(i),b(i));
    else
        b(i)=tan(pi/4-mod(y(i-1),2*pi));
        y(i)=phase1(2*y(i-1),b(i));  
    end
end

function y=phase1(w0,b0)
w=w0;b=b0;
y=w+2*atan(b*sin(w)./(1-b*cos(w)));

以上是实现附件里面的公式求出b值，另外还有一段英文说明。
这样应该得出的b值小于1才对，可是运行后的结果很怪，所以烦请高手指点一下程序问题所在，谢谢！