When I run this code in Matlab, it pops up "Function or variable 'six_bar' is not recognized”. How to solve the error?

clear;

i1=0.1;

i3=0.59;

i4=0.159;

i6=0.36;

i61=0.578;

omega1=0.87;

alpha=0;

hd=pi/180;

du=180/pi;

%调用子函数,six_bar返回牛头刨床的位移,角速度,角加速度

for n1 = 1:459

theta1(n1)=-2*pi+5.8119+(n1-1)*hd;

ii=[i1,i3,i4,i6,i61];

[theta,omega,alpha]=six_bar(theta1(n1),omega1,ii);

s3(n1)=theta(1);

theta3(n1)=theta(2);

theta4(n1)=theta(3);

sE(n1)=theta(4);

v2(n1)=omega(1);

omega3(n1)=omega(2);

omega4(n1)=omega(3);

vE(n1)=omega(4);

a2(n1)=alpha(1);

alpha3(n1)=alpha(2);

alpha4(n1)=alpha(3);

aE(n1)=alpha(4);

end

%3.位移,角速度，角加速度和牛头刨床图形输出

figure(1);

n1=1:459;

t=(n1-1)*2*pi/360;

subplot(2,2,1);%绘角位移及位移线图

plot(t,theta3*du,'r-.'); %'r-.'表示线条是红色的，并且线型是点画线

grid on;

hold on;

title('角位移及位移线图')

axis auto;

[haxes,~,~]=plotyy(t,theta4*du,t,sE); %

grid on;

hold on;

xlabel('时间/s');

axes(haxes(1));%haxes(1)表示左边的坐标轴的句柄,这行代码表示下面的代码是在左边的坐标轴上作图,修改

ylabel('角位移/\circ');

axes(haxes(2));

ylabel('位移/m');

hold on;

grid on;

text(1.15,-0.65,'\theta_3');

text(3.4,0.27,'\theta_4');

text(2.25,-0.15,'s_E');

subplot(2,2,2); %绘角速度及速度线图

plot(t,omega3,'r-.'); %'r-.'表示线条是红色的，并且线型是点画线

grid on;

hold on;

title('角速度及速度线图')

axis auto;

[haxes,~,~]=plotyy(t,omega4,t,vE); %

grid on;

hold on;

xlabel('时间/s');

axes(haxes(1));%haxes(1)表示左边的坐标轴的句柄,这行代码表示下面的代码是在左边的坐标轴上作图,修改

ylabel('角速度/rad\cdots∧{-1}');

axes(haxes(2));

ylabel('速度/m\cdots∧{-1}');

hold on;

grid on;

text(3.1,0.35,'\omega_3');

text(2.1,0.1,'\omega_4');

text(5.5,0.45,'v_E');

subplot(2,2,3); %绘角加速度及加速度线图

plot(t,alpha3,'r-.'); %'r-.'表示线条是红色的，并且线型是点画线

grid on;

hold on;

title('角加速度及加速度线图')

axis auto;

[haxes,~,~]=plotyy(t,alpha4,t,aE); %

grid on;

hold on;

xlabel('时间');

axes(haxes(1));%haxes(1)表示左边的坐标轴的句柄,这行代码表示下面的代码是在左边的坐标轴上作图,修改

ylabel('角加速度/rad\cdots∧{-2}');

axes(haxes(2));

ylabel('加速度/m\cdots∧{-2}');

hold on;

grid on;

text(1.5,0.3,'\alpha_3');

text(3.5,0.51,'\alpha_4');

text(1.5,-0.11,'a_E');

subplot(2,2,4); %牛头刨床机构

n1=20;

x(1)=0;y(1)=0;

x(2)=(s3(n1)*1000-50)*cos(theta3(n1));

y(2)=(s3(n1)*1000-50)*sin(theta3(n1));

x(3)=0;y(3)=i6*1000;

x(4)=i1*1000*cos(theta1(n1));

y(4)=s3(n1)*1000*sin(theta3(n1));

x(5)=(s3(n1)*1000+50)*cos(theta3(n1));

y(5)=(s3(n1)*1000+50)*sin(theta3(n1));

x(6)=i3*1000*cos(theta3(n1));

y(6)=i3*1000*sin(theta3(n1));

x(7)=i3*1000*cos(theta3(n1))+i4*1000*cos(theta4(n1));

y(7)=i3*1000*sin(theta3(n1))+i4*1000*sin(theta4(n1));

x(8)=i3*1000*cos(theta3(n1))+i4*1000*cos(theta4(n1))-900;

y(8)=i61*1000;

x(9)=i3*1000*cos(theta3(n1))+i4*1000*cos(theta4(n1))+600;

y(9)=i61*1000;

x(10)=(s3(n1)*1000-50)*cos(theta3(n1));

y(10)=(s3(n1)*1000-50)*sin(theta3(n1));

x(11)=x(10)+25*cos(pi/2-theta3(n1));

y(11)=y(10)-25*sin(pi/2-theta3(n1));

x(12)=x(11)+100*cos(theta3(n1));

y(12)=y(11)+100*sin(theta3(n1));

x(13)=x(12)-50*cos(pi/2-theta3(n1));

y(13)=y(12)+50*sin(pi/2-theta3(n1));

x(14)=x(10)-25*cos(pi/2-theta3(n1));

y(14)=y(10)+25*sin(pi/2-theta3(n1));

x(15)=x(10);

y(15)=y(10);

x(16)=0;

y(16)=0;

x(17)=0;

y(17)=i6*1000;

k=1:2;

plot(x(k),y(k));

hold on;

k=3:4;

plot(x(k),y(k));

hold on;

k=5:9;

plot(x(k),y(k));

hold on;

k=10:15;

plot(x(k),y(k));

hold on;

k=16:17;

plot(x(k),y(k));

hold on;

grid on;

axis ([-500 600 0 650]);

title('牛头刨床');

grid on;

xlabel('mm');

ylabel('mm');

plot(x(1),y(1),'o');

plot(x(3),y(3),'o');

plot(x(4),y(4),'o');

plot(x(6),y(6),'o');

plot(x(7),y(7),'o');

hold on;

grid on;

xlabel('mm');

ylabel('mm');

axis([-400 600 0 650]);

function [theta, omega, alpha]=six_bar(theta1,omega1,ii)

i1=ii(1);

i3=ii(2);

i4=ii(3);

i6=ii(4);

i61=ii(5);

%1.计算角位移和线位移

s3 =sqrt((i1*cos(theta1))*(i1*cos(theta1))+(i6+i1*sin(theta1))*(i6+i1*sin(theta1)));

%s3表示滑块2相对于CD杆的位移

theta3 =acos((i1*cos(theta1 ))/s3 );

%theta3表示杆3转过角度

theta4 =pi-asin((i61-i3*sin(theta3 ))/i4);

%theta4表示杆4转过角度

sE =i3*cos(theta3 )+i4*cos(theta4 );

%sE表示杆5的位移

theta(1)=s3;

theta(2)=theta3;

theta(3)=theta4;

theta(4)=sE;

%2.计算角速度和线速度

A=[sin(theta3 ),s3 *cos(theta3 ),0,0;

-cos(theta3 ),s3 *sin(theta3 ),0,0;

0,i3*sin(theta3 ),i4*sin(theta4 ),1;

0,i3*cos(theta3 ),i4*cos(theta4 ),0];

B=[i1*cos(theta1 );i1*sin(theta1 );0;0];

omega=A\(omega1*B);

v2 =omega(1);

%滑块2的速度　　

omega3 =omega(2);

%构件3的角速度　　

omega4 =omega(3);

%构件4的角速度

vE =omega(4);

%构件5的速度

%3.计算角加速度和加速度　　

A=[sin(theta3 ),s3 *cos(theta3 ),0,0;

%从动件位置参数矩阵

cos(theta3 ),-s3 *sin(theta3 ),0,0;

0,i3*sin(theta3 ),i4*sin(theta4 ),1;

0,i3*cos(theta3 ),i4*cos(theta4 ),0];

At=[omega3 *cos(theta3 ),(v2 *cos(theta3 )-s3*omega3 *sin(theta3 )),0,0;

-omega3 *sin(theta3 ),(-v2 *sin(theta3 )-s3 *omega3 *cos(theta3 )),0,0;

0,i3*omega3 *cos(theta3 ),i4*omega4 *cos(theta4 ),0;

0,-i3*omega3 *sin(theta3 ),-i4*omega4 *sin(theta4 ),0];

Bt=[-i1*omega1*sin(theta1 );-i1*omega1*cos(theta1 );0;0];

alpha=A\(-At*omega+omega1*Bt);

%机构从动件的加速度列阵　　

a2 =alpha(1);

%a2表示滑块2的加速度　　

alpha3 =alpha(2);

%alpha3表示杆3的角加速度

alpha4 =alpha(3);

%alpha4表示杆4的角加速度

aE =alpha(4);

%构件5的加速度

end

Answer 1

[Lone_wolf]

可以正常运行

clear;
i1=0.1;
i3=0.59;
i4=0.159;
i6=0.36;
i61=0.578;
omega1=0.87;
alpha=0;
hd=pi/180;
du=180/pi;
%调用子函数,six_bar返回牛头刨床的位移,角速度,角加速度
for n1 = 1:459
    theta1(n1)=-2*pi+5.8119+(n1-1)*hd;
    ii=[i1,i3,i4,i6,i61];
    [theta,omega,alpha]=six_bar(theta1(n1),omega1,ii);
    s3(n1)=theta(1);
    theta3(n1)=theta(2);
    theta4(n1)=theta(3);
    sE(n1)=theta(4);
    
    v2(n1)=omega(1);
    omega3(n1)=omega(2);
    omega4(n1)=omega(3);
    vE(n1)=omega(4);
    
    a2(n1)=alpha(1);
    alpha3(n1)=alpha(2);
    alpha4(n1)=alpha(3);
    aE(n1)=alpha(4);
end
    %3.位移,角速度，角加速度和牛头刨床图形输出
    figure(1);
    n1=1:459;
    t=(n1-1)*2*pi/360;
    subplot(2,2,1);%绘角位移及位移线图
    
    plot(t,theta3*du,'r-.');               %'r-.'表示线条是红色的，并且线型是点画线
    grid on;
    hold on;
    title('角位移及位移线图')
    axis auto;
    [haxes,~,~]=plotyy(t,theta4*du,t,sE); %
    grid on;
    hold on;
    
    xlabel('时间/s');
    axes(haxes(1));%haxes(1)表示左边的坐标轴的句柄,这行代码表示下面的代码是在左边的坐标轴上作图,修改
    ylabel('角位移/\circ');
    axes(haxes(2));
    ylabel('位移/m');
    hold on;
    grid on;
    text(1.15,-0.65,'\theta_3');
    text(3.4,0.27,'\theta_4');
    text(2.25,-0.15,'s_E');
    
    subplot(2,2,2);                         %绘角速度及速度线图
    plot(t,omega3,'r-.');               %'r-.'表示线条是红色的，并且线型是点画线
    grid on;
    hold on;
    title('角速度及速度线图')
    axis auto;
    [haxes,~,~]=plotyy(t,omega4,t,vE); %
    grid on;
    hold on;
    xlabel('时间/s');
    axes(haxes(1));%haxes(1)表示左边的坐标轴的句柄,这行代码表示下面的代码是在左边的坐标轴上作图,修改
    ylabel('角速度/rad\cdots∧{-1}');
    axes(haxes(2));
    ylabel('速度/m\cdots∧{-1}');
    hold on;
    grid on;
    text(3.1,0.35,'\omega_3');
    text(2.1,0.1,'\omega_4');
    text(5.5,0.45,'v_E');
    
    subplot(2,2,3);                         %绘角加速度及加速度线图
    plot(t,alpha3,'r-.');               %'r-.'表示线条是红色的，并且线型是点画线
    grid on;
    hold on;
    title('角加速度及加速度线图')
    axis auto;
    [haxes,~,~]=plotyy(t,alpha4,t,aE); %
    grid on;
    hold on;
    xlabel('时间');
    axes(haxes(1));%haxes(1)表示左边的坐标轴的句柄,这行代码表示下面的代码是在左边的坐标轴上作图,修改
    ylabel('角加速度/rad\cdots∧{-2}');
    axes(haxes(2));
    ylabel('加速度/m\cdots∧{-2}');
    hold on;
    grid on;
    text(1.5,0.3,'\alpha_3');
    text(3.5,0.51,'\alpha_4');
    text(1.5,-0.11,'a_E');
    
    
    subplot(2,2,4); %牛头刨床机构
    n1=20;
    x(1)=0;y(1)=0;
    x(2)=(s3(n1)*1000-50)*cos(theta3(n1));
    y(2)=(s3(n1)*1000-50)*sin(theta3(n1));
    x(3)=0;y(3)=i6*1000;
    x(4)=i1*1000*cos(theta1(n1));
    y(4)=s3(n1)*1000*sin(theta3(n1));
    x(5)=(s3(n1)*1000+50)*cos(theta3(n1));
    y(5)=(s3(n1)*1000+50)*sin(theta3(n1));
    x(6)=i3*1000*cos(theta3(n1));
    y(6)=i3*1000*sin(theta3(n1));
    x(7)=i3*1000*cos(theta3(n1))+i4*1000*cos(theta4(n1));
    y(7)=i3*1000*sin(theta3(n1))+i4*1000*sin(theta4(n1));
    x(8)=i3*1000*cos(theta3(n1))+i4*1000*cos(theta4(n1))-900;
    y(8)=i61*1000;
    x(9)=i3*1000*cos(theta3(n1))+i4*1000*cos(theta4(n1))+600;
    y(9)=i61*1000;
    x(10)=(s3(n1)*1000-50)*cos(theta3(n1));
    y(10)=(s3(n1)*1000-50)*sin(theta3(n1));
    x(11)=x(10)+25*cos(pi/2-theta3(n1));
    y(11)=y(10)-25*sin(pi/2-theta3(n1));
    x(12)=x(11)+100*cos(theta3(n1));
    y(12)=y(11)+100*sin(theta3(n1));
    x(13)=x(12)-50*cos(pi/2-theta3(n1));
    y(13)=y(12)+50*sin(pi/2-theta3(n1));
    x(14)=x(10)-25*cos(pi/2-theta3(n1));
    y(14)=y(10)+25*sin(pi/2-theta3(n1));
    x(15)=x(10);
    y(15)=y(10);
    x(16)=0;
    y(16)=0;
    x(17)=0;
    y(17)=i6*1000;
    k=1:2;
    plot(x(k),y(k));
    hold on;
    k=3:4;
    plot(x(k),y(k));
    hold on;
    k=5:9;
    plot(x(k),y(k));
    hold on;
    k=10:15;
    plot(x(k),y(k));
    hold on;
    k=16:17;
    plot(x(k),y(k));
    hold on;
    grid on;
    axis ([-500 600 0 650]);
    title('牛头刨床');
    grid on;
    xlabel('mm');
    ylabel('mm');
    plot(x(1),y(1),'o');
    plot(x(3),y(3),'o');
    plot(x(4),y(4),'o');
    plot(x(6),y(6),'o');
    plot(x(7),y(7),'o');
    hold on;
    grid on;
    xlabel('mm');
    ylabel('mm');
    axis([-400 600 0 650]);

    
function [theta, omega, alpha]=six_bar(theta1,omega1,ii)
    i1=ii(1);
    i3=ii(2);
    i4=ii(3);
    i6=ii(4);
    i61=ii(5);
    %1.计算角位移和线位移
    s3 =sqrt((i1*cos(theta1))*(i1*cos(theta1))+(i6+i1*sin(theta1))*(i6+i1*sin(theta1)));
    %s3表示滑块2相对于CD杆的位移
    theta3 =acos((i1*cos(theta1 ))/s3 );
    %theta3表示杆3转过角度
    theta4 =pi-asin((i61-i3*sin(theta3 ))/i4);
    %theta4表示杆4转过角度
    sE =i3*cos(theta3 )+i4*cos(theta4 );
    %sE表示杆5的位移
    theta(1)=s3;
    theta(2)=theta3;
    theta(3)=theta4;
    theta(4)=sE;
    %2.计算角速度和线速度
    A=[sin(theta3 ),s3 *cos(theta3 ),0,0;
    -cos(theta3 ),s3 *sin(theta3 ),0,0;
    0,i3*sin(theta3 ),i4*sin(theta4 ),1;
    0,i3*cos(theta3 ),i4*cos(theta4 ),0];
    B=[i1*cos(theta1 );i1*sin(theta1 );0;0];
    omega=A\(omega1*B);
    v2 =omega(1);
    %滑块2的速度　　
    omega3 =omega(2);
    %构件3的角速度　　
    omega4 =omega(3);
    %构件4的角速度
    vE =omega(4);
    %构件5的速度
    %3.计算角加速度和加速度　　
    A=[sin(theta3 ),s3 *cos(theta3 ),0,0;
    %从动件位置参数矩阵
    cos(theta3 ),-s3 *sin(theta3 ),0,0;
    0,i3*sin(theta3 ),i4*sin(theta4 ),1;
    0,i3*cos(theta3 ),i4*cos(theta4 ),0];
    At=[omega3 *cos(theta3 ),(v2 *cos(theta3 )-s3*omega3 *sin(theta3 )),0,0;
        -omega3 *sin(theta3 ),(-v2 *sin(theta3 )-s3 *omega3 *cos(theta3 )),0,0;
        0,i3*omega3 *cos(theta3 ),i4*omega4 *cos(theta4 ),0;
        0,-i3*omega3 *sin(theta3 ),-i4*omega4 *sin(theta4 ),0];
   Bt=[-i1*omega1*sin(theta1 );-i1*omega1*cos(theta1 );0;0];
   alpha=A\(-At*omega+omega1*Bt);
   %机构从动件的加速度列阵　　
   a2 =alpha(1);
   %a2表示滑块2的加速度　　
   alpha3 =alpha(2);
   %alpha3表示杆3的角加速度
   alpha4 =alpha(3);
    %alpha4表示杆4的角加速度
   aE =alpha(4);
    %构件5的加速度
end