我的OGY程序为啥没有控制住？

ieluxinhua

有谁做过OGY控制混沌程序的帮个忙看看我的程序为啥没有控制住，先谢谢了：

A=1.4;                              % 系统在A=1.4，B=0.3时呈现混沌运动
B=0.3;
d=0.08;
x(1)=0.32;                          % 系统初始值 x0=y0=0.32
xF=(B-1+((1-B).^2+4*A).^0.5)*0.5;   % 不动点（xF,xF)
f=[-2*xF B;1 0];                    % 系统在xF处的导数
g=((1-B).^2+4*A).^0.5*[1;1];        % 响应矢量矩阵
ju=-xF-(xF.^2+B).^0.5;              % 系统在xF点处不稳定方向上的特征根
eu=[1;-((xF.^2+B).^0.5-xF)/B];      % 系统在xF点处不稳定方向上单位向量
hu=[((xF.^2+B).^0.5+xF)/(2*((xF.^2+B).^0.5)) B/(-2*((xF.^2+B).^0.5))];%系统在不稳定方向上正交对偶矢量
for n=1:1000
    if n<2
        x(n+1)=x(n);
    else
        x(n+1)=A+B*x(n-1)-x(n).^2;      % 系统动力学方程   
        diag=[x(n)-xF;x(n)-xF];         % 控制的目的是使x(n+1)移到点xF的稳定方向上，即使diag在不稳定方向上的分量为零     
        xu=hu*diag;   
        P=(ju*hu*diag)/((ju-1)*(hu*g)); % 控制率表达
        if -(1-1/ju)*hu*g*d<xu<(1-1/ju)*hu*g*d  % 控制启动条件
           x(n+1)=A+B*x(n-1)-x(n).^2+P;
        else
           x(n+1)=A+B*x(n-1)-x(n).^2;
        end
     end
end
n=1:1000;
plot(n,x(n),'r.');
grid on;