[leek008]

主程序
y0=[0;0;0;0];t0=0;tf=3;t=[t0,tf];
while t0<=tf
if y0(1)<=0
[t,y]=ode45('contacty',t,y0);
n=length(y(:,1));
for i=1:n;
B=y(i,1);
if B>0
break
end
end
plot(t(1:i),y(1:i,1))
hold on
t0=t(i);y0=[y(i,1);y(i,2);y(i,3);y(i,4)];
else
[t,y]=ode45('jumpy',t,y0);
n=length(y(:,1));
for i=1:n;
B=y(i,1);
if B<=0
break
end
end
plot(t(1:i),y(1:i,1))
hold on
t0=t(i),y0=[y(i,1);y(i,2);y(i,3);y(i,4)];
end
end
m函数
function dydt=contacty(t,y)
dydt=[y(2);(20000*sin(40*t)-200*y(2)-4040000*y(1)+100*y(4)+3200000*y(3))/100;y(4);(-100*y(4)-3200000*y(3)+100*y(2)+3200000*y(1))/100];

function dydt=jumpy(t,y)
dydt=[y(2);(2000*sin(40*t)-1000-100*y(2)-840000*y(1)+100*y(4)+840000*y(3))/100;y(4);(100*y(2)+840000*y(1)-100*y(4)-840000*y(3))/100];