a=0;
b=pi/2;
a1=a;
a3=b;
a2=a;
a4=b;
while a3-a1>0.01
f1=-atan(2*sin(2*a1)/(5-cos(2*a1)));
f3=-atan(2*sin(2*a3)/(5-cos(2*a3)));
a2=(a1+a3)/2;
f2=-atan(2*sin(2*a2)/(5-cos(2*a2)));
c1=(f3-f1)/(a3-a1);
c2=((f2-f1)/(a2-a1)-c1)/(a2-a3);
if c2~=0
a4=0.5*(a1+a3-c1/c2);
f4=-atan(2*sin(2*a4)/(5-cos(2*a4)));
if (4-a1)*(a3-a4)>0
    if a2<a4
        if f2<f4
            a3=a4;  
        else
            a1=a2;
        end
    else
        if f2<f4
            a1=a4;
        else
            a3=a2;
        end
    end
end
end
end
if f2<f4
    x=a2;
else
    x=a4;
end
f=-atan(2*sin(2*x)/(5-cos(2*x)));
x=x*180/pi;
f=f*180/pi;
fprintf('最大壓力角α=%d/n,相應(yīng)的凸輪轉(zhuǎn)角φ=%d/n',f,x)