matlab画频散曲线,关于lamb频散曲线的绘制问题

clc

clear

close all

n1=0;

syms c;   %定义速度c为符号函数

err=[];      %xx，err存贮误差

eps=1e-12;         %精度

ct=3.130;cl=6.350;%铝板的横波速度3130m/s，纵波速度6350m/s   这边的速度最好是给出的个位数，如果太大，完蛋！那tanh会趋近于0，所以图形就不对

d=1;

m=0.5:0.1:3.1;

A=[];

cp1 = [];

n1 = 0;

k = 1;

yy = [];

for j=1:length(m)

xx=[];     %%存贮当前cp下频率f的解

cp=m(j);

f1=0.1:0.01:20;

f=0.1;

y1=((cp/ct)^2-2)^2*tanh((d*pi*f/cp)*sqrt(1-(cp/cl)^2));

y2=4*sqrt(1-(cp/ct)^2)*sqrt(1-(cp/cl)^2)*tanh((d*pi*f/cp)*sqrt(1-(cp/ct)^2));

%plot(f,y1,f,y2)

b1=subs(y1-y2);

y11=[y1];

y22=[y2];

for i=2:length(f1)

f=f1(i);

y1=((cp/ct)^2-2)^2*tanh((d*pi*f/cp)*sqrt(1-(cp/cl)^2));

y2=4*sqrt(1-(cp/ct)^2)*sqrt(1-(cp/cl)^2)*tanh((d*pi*f/cp)*sqrt(1-(cp/ct)^2));

y11=[y11;y1];

y22=[y22;y2];

b2=subs(y1-y2);

if  b2==0

xx=[xx;f];

else if b1*b2<0     %在b1*b2<0是可以改进，取(f1+f2)/2

%  xx=[xx;(f1(i)+f1(i-1))/2];

a = f1(i-1);b = f1(i);

while ((b-a)>eps)

c = (b+a)/2;

fa = ((cp/ct)^2-2)^2*tanh((d*pi*a/cp)*sqrt(1-(cp/cl)^2))-4*sqrt(1-(cp/ct)^2)*sqrt(1-(cp/cl)^2)*tanh((d*pi*a/cp)*sqrt(1-(cp/ct)^2));

fb = ((cp/ct)^2-2)^2*tanh((d*pi*b/cp)*sqrt(1-(cp/cl)^2))-4*sqrt(1-(cp/ct)^2)*sqrt(1-(cp/cl)^2)*tanh((d*pi*b/cp)*sqrt(1-(cp/ct)^2));

fc = ((cp/ct)^2-2)^2*tanh((d*pi*c/cp)*sqrt(1-(cp/cl)^2))-4*sqrt(1-(cp/ct)^2)*sqrt(1-(cp/cl)^2)*tanh((d*pi*c/cp)*sqrt(1-(cp/ct)^2));

if fa*fc < 0

b=c;

else

a=c;

end

end

xx = [xx;c];

end

end

b1=b2;

end

n=length(xx);

%if n>n1

%   n1=n;