求matlab高手指导程序错误

Question

function four_bar
global L0 L1 L2 L3 th1;
%str=inputdlg({'please input rack_length L0:','please input crank_length L1:',
 %  'please input linkage_length L2:','please input rock_length L3:',
  % 'please input w1(rad/sec):','please input repeat times m:',
   %'please input initial theta3:'},'please input!',1,{'20','8','25','20','100','10','1'});
%clc;
L0=20;
L1=8;
L2=25;
L3=20;
w1=100;
m=10;
theta3=1;
t=linspace(0,2*pi/w1,181);
c = [t(:,1:1-1) t(:,1+1:end)];
b = [c(:,1:1-1) c(:,1+1:end)];
dt=2*pi/w1/180;
theta1=w1*t;
th1=theta1(1);
%y=four_bar_equation(x);
theta3(1)=fzero(@four_bar_equation,theta3);
x_a=0;y_a=0;
x_b=x_a+L1*cos(theta1(1));y_b=y_a+L1*sin(theta1(1));
x_d=x_a+L0;y_d=y_a;
x_c=x_d+L3*cos(theta3(1));y_c=y_d+L3*sin(theta3(1));
figure('name','four_bar_move...','numbertitle','off','position',[232 199 617 479]);
line([x_a;x_d],[y_a;y_d],'color','m','LineStyle','-','linewidth',3,'erasemode','none');
line([x_a;x_a],[y_a;y_a],'color','black','LineStyle','.','erasemode','xor','MarkerSize',30);
line([x_d;x_d],[y_d;y_d],'color','black','linestyle','.','erasemode','xor','MarkerSize',30);
dot_b=line([x_b;x_b],[y_b;y_b],'color','black','linestyle','.','erasemode','xor','MarkerSize',25);
dot_c=line([x_c;x_c],[y_c;y_c],'color','black','linestyle','.','erasemode','xor','MarkerSize',25);
line_ab=line([x_a;x_b],[y_a;y_b],'color','r','linestyle','-','linewidth',3,'erasemode','xor');
line_bc=line([x_b;x_c],[y_b;y_c],'color','g','linestyle','-','linewidth',3,'erasemode','xor');
line_cd=line([x_c;x_d],[y_c;y_d],'color','b','linestyle','-','linewidth',3,'erasemode','xor');
axis([-10,30,-10,25]);
%axis equal;
for i=1:m
    for j=2:181
        th1=theta1(j);
        theta3(j)=fzero(@four_bar_equation,theta3(j-1));
        x_b=x_a+L1*cos(th1);y_b=y_a+L1*sin(th1);
        x_b_c1(j,1)=x_b;y_b_c1(j,1)=y_b;
        x_c=x_d+L3*cos(theta3(j));
        y_c=y_d+L3*sin(theta3(j));
        x_c_c1(j,1)=x_c;y_c_c1(j,1)=y_c;
        set(dot_b,'xdata',x_b,'ydata',y_b);
        set(dot_c,'xdata',x_c,'ydata',y_c);
        set(line_ab,'xdata',[x_a;x_b],'ydata',[y_a;y_b]);
        set(line_bc,'xdata',[x_b;x_c],'ydata',[y_b;y_c]);
        set(line_cd,'xdata',[x_c;x_d],'ydata',[y_c;y_d]);
        drawnow;
    end
end
w3=diff(theta3)/dt;
a3=diff(w3)/dt;
subplot(1,3,1);
plot(t,theta3);
subplot(1,3,2);
plot(c,w3);
subplot(1,3,3);
plot(b,a3);
function y=four_bar_equation(x)
global L0 L1 L2 L3 th1;
%y=2*x-1;
y=L1*cos(th1)+L2*sqrt(1-(L3*sin(x)-L1*sin(th1))^2/L2/L2)-L3*cos(x)-L0;