我想画一个刚度曲线图，但是我编的程序有问题，我也不知道那错了。就是图画不出来。帮我看看吧，谢谢！

这是我编写的程序

z1 = 25;                  %齿数

z2 = 35;

m = 4;                    %模数

a = deg2rad(20);     %齿轮压力角

c = 0.25;                 %顶隙系数

ha = 1;                   %齿顶高系数

x = 0;                     %变位系数

E = 2*10^9;           %弹性模量

v = 0.3;                  %泊松比

L = 40;                   %齿宽

rb1 = m*z1*cos(a);

rb2 = m*z2*cos(a);           %基圆半径

ra1 = m*(z1+2*ha);

ra2 = m*(z2+2*ha);                  %齿顶圆半经

inva = tan(a)-a;                         %渐开线函数

S = (pi/2+2*x*tan(a))*m;          %分度圆齿厚

Sb1 = cos(a)*(S+m*z1*inva);    %基圆齿厚

Sb2 = cos(a)*(S+m*z2*inva);

kh = pi*E*L/(4*(1-v^2));          %赫兹接触刚度

alphak2 = Sb1/(2*rb1);

rho1 = 34.2020;

alphak1 = rho1/rb1-alphak2;

alphak4 = Sb2/(2*rb2);

rho2 = 47.8828;

alphak3 = rho2/rb2-alphak4;

syms alpha1 beta1 gamma1 alpha2 beta2 gamma2;

f1 = 3*(1+cos(alphak1)*((alphak2-alpha1)*sin(alpha1)-cos(alpha1)))^2*(alphak2-alpha1)*cos(alpha1) ...

/(2*E*L*(sin(alpha1)+(alphak2-alpha1)*cos(alpha1)))^3;

f2 = 1.2*(1+v)*(alphak2-beta1)*cos(beta1)*(cos(alphak1))^2*alphak1 ...

/(E*L*(sin(beta1)+(alphak2-beta1)*cos(beta1)));

f3 = (alphak2-gamma1)*cos(gamma1)*(sin(alphak1))^2/(2*E*L*(sin(gamma1)+(alphak2-gamma1)*cos(gamma1)));

f4 = 3*(1+cos(alphak3)*((alphak4-alpha2)*sin(alpha2)-cos(alpha2)))^2*(alphak4-alpha2)*cos(alpha2) ...

/(2*E*L*(sin(alpha2)+(alphak4-alpha2)*cos(alpha2)))^3;

f5 = 1.2*(1+v)*(alphak4-beta2)*cos(beta2)*(cos(alphak3))^2*alphak3 ...

/(E*L*(sin(beta2)+(alphak4-beta2)*cos(beta2)));

f6 = (alphak4-gamma2)*cos(gamma2)*(sin(alphak3))^2/(2*E*L*(sin(gamma2)+(alphak4-gamma2)*cos(gamma2)));

inv1 = int(f1,alpha1,-alphak1,alphak2);

inv2 = int(f2,beta1,-alphak1,alphak2);

inv3 = int(f3,gamma1,-alphak1,alphak2);

inv4 = int(f1,alpha2,-alphak3,alphak4);

inv5 = int(f2,beta2,-alphak3,alphak4);

inv6 = int(f3,gamma2,-alphak3,alphak4);

kb1 = 1/inv1; %弯曲刚度

ks1 = 1/inv2; %剪切刚度

kc1 = 1/inv3; %轴向压缩刚度

kb2 = 1/inv4; %弯曲刚度

ks2 = 1/inv5; %剪切刚度

kc2 = 1/inv6; %轴向压缩刚度

K = 1/(1/kh+1/kb1+1/kb2+1/ks1+1/ks2+1/kc1+1/kc2)；

plot(K)

运行结果为：

错误使用 plot

数据必须为可转换为双精度值的数值、日期时间、持续时间或数组。

出错 Kt (line 62)

plot(K)

2018-4-5 21:08 上传

只是主要的公式