请教：频率特性曲线绘制程序代码有什么问题？

baipi

附件是特性关系方程
程序代码为
clear;
omega=100;%扰频率
p=20;%外力幅值
h=0.1e-3;%板厚
E=70e9;%杨氏摸量
ro=2700;%密度
v=0.33;%泊松比
sigma=3.6e7;%电导率
BZ=.1;%磁场强度
b=pi/8;
a=pi/6;
DN=E*h/(1-v^2);%抗拉刚度
DM=E*h^3/(12*(1-v^2));%抗弯刚度
d1=ro*h;%加速度前系数
a1=(pi/a)^2+(pi/b)^2;
a2=9*a1^2-16*(pi/a)^2*(pi/b)^2;
d3=a1^2*DM;
d4=a2*DN/32;
beta1=d4/d1*h^2;%三次系数
alfa1=d3/d1;%固有频率平方相
omga=omega/sqrt(alfa1);
    f1=p;%力幅
   f2=f1/d1;
   f=f2*(sqrt(abs(beta1/alfa1))/alfa1);%外力 *
  d21=sigma*h^3*BZ^2/12;
d2=d21*a1;%速度前系数
mu=d2/d1*sqrt(alfa1);%阻尼系数*sqrt(beta1)/alfa1
B=f;A1=[];z1=[];S=[];AA=[];
i=0;j=0;
for A=0.01:0.001:10;
    z=mu^4/4-mu^2*(1-3*beta1*A.^2/4)+(B./A).^2;
    if z>=0
        i=i+1;
        A1(i)=A;
        z1(i)=z;
    end
end
A=A1;
z=z1;
s1=(sqrt(1-mu^2/2-3*beta1*A.^2/4+sqrt(z1)));
ceta=atan(-mu*s1./(1-s1.^2-3*beta1*A1.^2/4));
ceta1=mod(ceta,pi);
%plot(s1,A)
i=0;
for A2=A;
    i=i+1;
    ss=1-mu^2/2-3*beta1*A.^2/4-sqrt(z(i));
    if ss>=0
        j=j+1;
        S(j)=sqrt(ss);
        AA(j)=A2;
    end
end
ceta=atan(-mu*S./(1-S.^2+3*beta1*AA.^2/4));
ceta=mod(ceta,pi);
hold on;
figure (1)
plot(S,AA,'k',s1,A,'k')          %幅频
figure (2)
plot(S,ceta,'k',s1,ceta1,'k');   %相频

baipi

怎么没有人指教啊？？？？

baipi

问题已经解决了

十全武功

怎么解决的啊？