function [y1,y2,y3]=OdePP()    %画函数图像
PeriodT=0.5;p=2; %脉冲周期、脉冲量
NumOfStep=50;NumOfPer=4; %NumOfStep是每半个周期取多少点画图,NumOfPer 是画多少个周期的相图
Init_y1=1*10^6;Init_y2=1*10^3;Init_y3=0;Init_t=0;  %初值和初始时刻
NStep1=fix(NumOfStep*1);
NStep2=NumOfStep-NStep1;
y1=zeros((NumOfStep+2)*NumOfPer,1);%存放第一维数据的数组
y2=zeros((NumOfStep+2)*NumOfPer,1);%存放第二维数据的数组
y3=zeros((NumOfStep+2)*NumOfPer,1);%存放第三维数据的数组
t=zeros((NumOfStep+2)*NumOfPer,1);%存放时间值数据的数组
ppp=odeset('RelTol',1e-1,'AbsTol',1e-1);  %控制误差的参数
for i=1:NumOfPer     %循环一次画一个周期的轨线
     if i==1
       y1(1)=Init-y1;y2(1)=Init-y2;y3(1)=Init-y3 ;
     else   
         %下面3个语句确定在T时刻的脉冲值
       Init-y1=y1*((i-1)*(NumOfStep+2));
       Init-y2=y2*((i-1)*(NumOfStep+2));  
       Init-y3=y3*((i-1)*(NumOfStep+2))*(1-p);   
       y1*((i-1)*(NumOfStep+2)+1)=Init-y1;
       y2*((i-1)*(NumOfStep+2)+1)=Init-y2;
       y3*((i-1)*(NumOfStep+2)+1)=Init-y3;
end
     t((i-1)*(NumOfStep+2)+1)=(i-1)*PeriodT+Init-t;
     sol=ode45(@OdePP,[(i-1)*PeriodT+Init-t,(i-1+1)*PeriodT+Init-t],[Init-y1,Init-y2,Init-y3],ppp);
     j=linspace(PeriodT*(i-1)+Init-t,(i-1+1)*PeriodT+Init-t,NStep1);
     y1(((i-1)*(NumOfStep+2)+2):((i-1)*(NumOfStep+2)+1+NStep1))=deval(sol,j,1);
     y2(((i-1)*(NumOfStep+2)+2):((i-1)*(NumOfStep+2)+1+NStep1))=deval(sol,j,2);
     y3(((i-1)*(NumOfStep+2)+2):((i-1)*(NumOfStep+2)+1+NStep1))=deval(sol,j,3);
     t(((i-1)*(NumOfStep+2)+2):((i-1)*(NumOfStep+2)+1+NStep1))=j;
        %以上第一次脉冲完成   
end
save TrajData  %此语句用于保存所有内存数据
figure(1);     %此语句绘出相图
  plot3(y1,y2,y3,'k');
   box on   xlabel('y-1');ylabel('y-2');zlabel('y-3');
   hgsave('Traj');%此语句用于保存绘出的图
  figure(2);
   plot(t,y1,'k');
   xlabel('t');ylabel('y-1') ;
    hgsave('tseries-y1');%绘制并保存第一维的时间序列图
 function dydt=OdePP(t,y)  %决定方程的函数
r=1; k=3;a=3;b=4; d=0.2;p=0.4 ; q=0.45; e=0.2;  %系统参数
dy = zeros(3,1);
dy(1) =r*y(1)*(1-y(1)/k)-a*y(1)*y(2);
dy(2) =y(2)*((b-1)*a*y(1)-d-p*y(3));
dy(3) =y(3)*(q*y(2)-e);
   请问谁做过脉冲微分方程的数值模拟？以上是我根据网上的资源编的脉冲程序，不过有问题运行不出来，请有经验的高手帮我看看，这里先谢谢啦！如果成功的话，有重酬！

[此贴子已经被作者于2018-7-14 17:41编辑过]