function newx=Henon(x);
% Henon方程[差分方程]
% newx=Henon([x;y;p;q])
%
% 方程如下:
%   x(k+1)=-p*x(k)*x(k)+y(k)+1
%   y(k+1)=q*x(k)%% Example(函數圖象):
%
%   X=[1;0;1.4;.3];Y=[];
%   for i=1:1000
%       X=feval(@Henon,X);
%       Y(i,:)=X(1:2,1);
%   end
%   plot(Y(:,1),Y(:,2),'.');
%   title('Henon映射圖'),xlabel('x'),ylabel('y')
%
% Example(分岔圖):
% 
%   Z=[]; 
%   for p=linspace(0,1.4,300);
%       x=[1;0;p;.3]; 
%       for k=1:300; 
%           x=Henon(x);
%           if k>60
%             Z=[Z,p+x(1)*i]; 
%           end
%       end 
%   end 
%   plot(Z,'.','markersize',1) 
%   title('Henon映射分岔圖'),xlabel('p'),ylabel('x')
%
% Example(最大Lyapunov指數譜圖):
%
%   d0=1e-8;
   Z=[]; 
%   for p=linspace(0,1.4,300)
%       le=0;
%       lsum=0;
%       x=[1;0;p;.3];
%       x1=[1;d0;p;.3];
%       for k=1:800
%           x=Henon(x);
%           x1=Henon(x1);
%           d1=sqrt((x(1)-x1(1))^2+(x(2)-x1(2))^2);
%           x1=x+(d0/d1)*(x1-x);
%           if k>100
%               lsum=lsum+log(d1/d0);
%           end
%       end
%       le=lsum/(k-100);
%       Z=[Z,p+le*i]; 
%   end
%   plot(Z,'-') 
%   title('Henon最大Lyapunov指數譜圖'),xlabel('p'),ylabel)('lyapunov')
%   grid on
%
% Author：yujunjie
% Author's email: ustb03-07@yahoo.com.cn
newx(1,1)=-x(3)*x(1)*x(1)+x(2)+1;
newx(2,1)=x(4)*x(1);
newx(3,1)=x(3);
newx(4,1)=x(4);