Hopf分岔点搜索、同宿分岔，高手帮忙看看，是不是程序哪里有问题？

gouxu

本来是想做一个搜索Hopf分岔点的程序，结果在分岔参数Q1=11.14时，状态方程的解出现NAN的情况，和实例对比，发现在Q1=11.2的时候发生同宿分岔。百思不得其解。对照了实例，没有发现有输入错误的地方。特求助各位达人！望给予指导！
用matcont做了仿真，发现Q1=11.14没有什么异常情况，和实例吻合的很好。
实例中，(X,Q1)=[0.3090,0,0.12,1.1,10.945] 为Hopf点
              (X,Q1)=[0.3465,0,0.1374,0.9294,11.405] 为Hopf点
              (X,Q1)=[0.3485,0,0.1384,0.9208,11.411] 为SNB点

///////////////////////////////////////////////////////非线性微分方程组
function F=myfun(t,x,Q1)

F=[x(2);
   16.6667*sin(x(3)-x(1)+0.0873)*x(4)-0.1667*x(2)+1.8807;
   496.8718*x(4)^2-166.6667*cos(x(3)-x(1)-0.0873)*x(4)-666.6667*cos(x(3)-0.2094)*x(4)-93.3333*x(4)+33.3333*Q1+43.3333;
   -78.7638*x(4)^2+26.2172*cos(x(3)-x(1)-0.0124)*x(4)+104.8689*cos(x(3)-0.1346)*x(4)+14.5229*x(4)-5.2288*Q1-7.0327];

///////////////////////////////////////////////////一下是主程序
%搜索Hopf分岔点
tic
A0=eye(4);
A=j*100*pi*A0;
h=1e-2;
t0=0;   %开始时间
tf=15;  %终了时间
x0=[0.3;0;0.2;0.95]; %初值
Q1=10:0.01:12;  %分岔参数取值范围
for k=1:length(Q1) % length(Q1)为Q1的点数
    kk=Q1(k); %每一个Q1值赋给kk
    for t=t0:h:tf   %四阶龙格库塔法解方程组
        k1=myfun(t,x0,kk);
        k2=myfun(t+h/2,x0+h/2*k1,kk);
        k3=myfun(t+h/2,x0+h/2*k2,kk);
        k4=myfun(t+h,x0+h*k3,kk);
        x=x0+h/6*(k1+2*k2+2*k3+k4);
        x0=x;
    end
    x1(k)=x(1);
    x2(k)=x(2);  %把每一个Q1值下的解保存
    x3(k)=x(3);
    x4(k)=x(4);
   
    J=[0,1,0,0;
       -16.6667*cos(x(3)-x(1)+0.0873)*x(4),-0.1667,16.6667*cos(x(3)-x(1)+0.0873)*x(4),16.6667*sin(x(3)-x(1)+0.0873);
       -166.6667*sin(x(3)-x(1)-0.0873)*x(4),0,166.6667*sin(x(3)-x(1)-0.0873)*x(4)+666.6667*sin(x(3)-0.2094)*x(4),496.8718*2*x(4)-166.6667*cos(x(3)-x(1)-0.0873)-666.6667*cos(x(3)-0.2094)-93.3333;
       26.2172*sin(x(3)-x(1)-0.0124)*x(4),0,-26.2172*sin(x(3)-x(1)-0.0124)*x(4)-104.8689*sin(x(3)-0.1346)*x(4),-78.7638*2*x(4)+26.2172*cos(x(3)-x(1)-0.0124)+104.8689*cos(x(3)-0.1346)+14.5229];
     
    %构造In,m
    if det(A-J)<=1e-3  && myfun(t,x,Q1)<=1e-3    %判断是否为Hopf分岔点
        XHopf=x;
        disp('Q1=')
        disp(Q1(k))
        disp('Hopf Bifurcation XHopf=')
        disp(XHopf)
    end
   
end
plot(Q1,x1,'-b',Q1,x4,'-r',Q1,x3,'-y');
toc

////////////////////////////////////////////////////////////////////
Hopf分岔点判据：(见《基于分岔理论的电力系统电压稳定性分析》一书P12)
F(X,Q1)=0 && det(j*2*pi*f*A)-J(X,Q1))=0
J为雅克比矩阵

上式中，A=eye(4);
//////////////////////////////////////////////////////////////////
附图：
红色为x(4)   蓝色为x(1)   黄色为x(3)

gouxu

这是用matcont仿真出来的图形：

求救人

你是用A MODEL OF VOLTAGE COLLAPSE IN ELECTRIC POWER SYSTEMS它做的吧。