有偿改MATLAB代码

上图为预期目标，可以不要简图

现有代码：

clear;
l1 = 100;  % 主动杆杆长
l2 = 100;  % 连杆杆长
hd = pi/180;
du = 180/pi;
theta1=30*du;
theta2=30*du;

omega1 = 30; % 主动杆角速度
alpha1 = 0; % 主动杆角加速度

s=l2 * cos(theta2) -l1 * cos(theta1); 
A=[l2*sin(theta2), 1; l2*cos(theta2), 0]; 
B=[-l1*sin(theta1); l1*cos(theta1)];  
% ‘\’求解关于 x 的线性方程组 Ax = B 
%  eg: x=A\B 
omega=A \ (omega1*B);         
omega2=omega(1); 
v=omega(2);

Aa=[omega2*l2*cos(theta2), 0; 
    -omega2*l2*sin(theta2), 0];        
Ba= [omega1*l1*cos(theta1);       
    -omega1*l1*sin(theta1)];      
alpha=A \(-Aa*omega+alpha1*B+omega1*Ba);  
alpha2 = alpha(1); 
a = alpha(2); 

% for循环获取0-4pi构件2角位移、滑块位移、构件2角速度、滑块线速度、构件2角加速度、滑块线加速度
for n1= 1:720
    theta1(n1)=(n1-1)*hd;  % 0-4pi
    [theta2(n1) , s(n1) , omega2(n1), v(n1) , alpha2(n1), a(n1)] = slider_crank(theta1(n1) , omega1, alpha1, l1, l2);
end
figure(1);
n = 1:720;
% 位移曲线图
subplot(2, 2, 1);
yyaxis left;
plot(theta1*du, theta2*du);
ylabel('连杆角位移/(\circ)');
yyaxis right;
plot(theta1*du, s);
ylabel('滑块位移/mm');
xlabel('曲柄转角\theta_1/(\circ))');
title('位移曲线图');
xlim([0, 720]);
x_new = linspace(0, 720, 5);
xticks(x_new)
grid on;

%速度曲线图
subplot(2,2,2);
yyaxis left;
plot(theta1*du, omega2);
ylabel('连杆角速度/(rad/s)');
yyaxis right;
plot(theta1*du, v/1000);
ylabel('滑块速度/(m/s)');
xlabel('曲柄转角\theta_1/(\circ))');
xlim([0, 720]);
x_new = linspace(0, 720, 5);
xticks(x_new)
title('速度曲线图');
grid on;
%加速度曲线图
subplot(2,2,3);
yyaxis left;
plot(theta1*du, alpha2);
ylabel('连杆角加速度/rad\cdots^{-2}');
yyaxis right;
plot(theta1*du, a/10^6);
ylabel('滑块加速度/m\cdots^{-1}');
xlabel('曲柄转角\theta_1/(\circ))');
xlim([0, 720]);
x_new = linspace(0, 720, 5);
xticks(x_new)
title('加速度曲线图');
grid on;
function [theta2,s,omega2,v,alpha2,a]=slider_crank(theta1,omega1,alpha1,l1,l2)
sa=l2 * cos(theta2) - l1 * cos(theta1);
A=[l2*sin(theta2), 1; l2*cos(theta2), 0]; 
B=[l1*sin(theta1); l1*cos(theta1)];  
% ‘\’求解关于 x 的线性方程组 Ax = B 
%  eg: x=A\B 
omega=A \ (omega1*B);         
omega2=omega(1); 
v=omega(2); 
Aa=[omega2*l2*cos(theta2), 0; 
    -omega2*l2*sin(theta2), 0];        
Ba= [omega1*l1*cos(theta1);       
    -omega1*l1*sin(theta1)];      
alpha=A \(-Aa*omega+alpha1*B+omega1*Ba);  
alpha2 = alpha(1); 
a = alpha(2); 
end
现存在问题：