被动行走模型动力学分析的matlab实现

%% 定义计算的步长, 设置变量的初始值
clear;clc;close all
Delta=0.001;     % 步长
t=0:Delta:8;     % 定义自变量 t
n=length(t); 
Y(:,1)=[0.2;-0.2;-1.5-1.0];  % 变量 th1 th2 th1d th2d 的初始值

%% 自定义龙格库塔法, 求解微分方程组
for k=1:n-1
    z1=f(t(k),Y(:,k));
    z2=f(t(k)+Delta/2,Y(:,k)+z1*Delta/2);
    z3=f(t(k)+Delta/2,Y(:,k)+z2*Delta/2);
    z4=f(t(k)+Delta,Y(:,k)+z3*Delta);
    Y(:,k+1)=Y(:,k)+Delta*(z1+2*z2+2*z3+z4)/6;
end
th1=Y(1,:);
th2=Y(2,:);
th1d=Y(3,:);
th2d=Y(4,:);

%% 绘制 x y 的求解结果
figure
set(gcf,'units','normalized','position',[0.15 0.2 0.7 0.6]);  % 设置 figure 窗口的位置和尺寸
subplot(1,2,1)
plot(t,th1,'b',th2,'r')
xlabel('t')
ylabel('th1','th2')

subplot(1,2,2)
plot(t,th1d,'b',th2d,'r')
xlabel('t')
ylabel('th2d','th2d')
  
索引超出数组元素的数目(3)。

出错 f (line 3)
th1=Y(1);th2=Y(2);th1d=Y(3);th2d=Y(4);

出错 RungeKutta (line 10)
    z1=f(t(k),Y(:,k));

function F=f(t,Y)
m1=1;a1=2;m2=1;a2=1;g=9.81;c=2;
th1=Y(1);th2=Y(2);th1d=Y(3);th2d=Y(4);
f1=m1*(a1*a2-a1*a1)*sin(th1+th2)*th2d*(th2d-2*th1d)-m1*g*(a1+a2)*sth1*cos(c)-m2*g*a1*cos(c)*sth1;
f2=m1*(a1*a2-a1*a1)*sin(th1+th2)*th1d*th1d+m1*g*(a1-a2)*sth2*cos(c);
F=[f1;f2];
end

双足机器人二维被动行走模型的构建，运动微分方程已定义，运用matlab利用runge-kutta法跑出角度、角速度与时间t的关系图像。

你好，我是有问必答小助手。为了技术专家团更好地为您解答问题，烦请您补充下（1）问题背景详情，（2）您想解决的具体问题，（3）问题相关代码图片或者报错信息。便于技术专家团更好地理解问题，并给出解决方案。

您可以点击问题下方的【编辑】，进行补充修改问题。