[求助]poincare截面程序求助

[zhailiangjun]

请各位大虾帮忙看一下我写的这个计算poincare截面图的程序究竟出了什么问题呐。就是画不出图来，自己查了很久了，可还是找不到错误，也不知道该怎样改，小弟先谢谢了。

function dxdt=zhengzefun(t,x)
global w;
global f;
global j;
w=1;
f=0.1;
j=4.5;
dxdt=zeros(4,1);
%x(1)=p1 x(2)=q1 x(3)=p2 x(4)=q2
dxdt(1)=-w*x(2)-2*f*((4*j-x(2)^2-x(1)^2)/4/j)^(1/2)*x(4)-f*x(2)^2*x(4)*((4*j-x(2)^2-x(1)^2)/4/j)^(-1/2)/j;
dxdt(2)=w*x(1)-x(1)*x(2)*x(4)*((4*j-x(2)^2-x(1)^2)/4/j)^(-1/2)/j;
dxdt(3)=-w*x(4)-2*f*((4*j-x(2)^2-x(1)^2)/4/j)^(1/2)*x(2);
dxdt(4)=w*x(3);

上面是一个哈密顿正则运动方程，我试图画出它的poincare截面图。所以参考本论坛上许多类似的帖子，然后编出了下面的程序。

%poincare
%x(1)=p1 x(2)=q1 x(3)=p2 x(4)=q2 lamt=f
global w,
global f,
global j
w=1;
f=0.5;
j=4.5;
[t,x]=ode113('zhengzefun',[0:6.3021e-004:630.211],[2.79429;-0.4326;0.1253;-1.6656]);
for i=1:1length(x)-1)
if x(i,4)<0
if x(i+1,4)>0;
if x(i,3)>0
qq(i)=x(i,1);
pp(i)=x(i,2);
end
end
end
if x(i,4)>0
if x(i+1,4)<0;
if x(i,3)>0
qq(i)=x(i,1);
pp(i)=x(i,2);
end
end
end
if x(i,4)=0
if x(i,3)>0
qq(i)=x(i,1);
pp(i)=x(i,2);
end
end
end
plot(qq(,pp(,'.')

我是想画出x（4）=0的面上，当x（3）>0时的x（1），x（2）的截面图，（x（1）与x（2）的数值校正没有写出），但是每次都得出了很少的点。改变了积分步长和积分时间还是不行。是不是我对程序理解的不对呐，请教大家了，多谢多谢。