function henonphase
 mysiz=500;
 clf
 e=1/12;
 hold on
 Himage = ones(mysiz,mysiz,3);
 p=.5; q=0; r=.05; s=100;
 Himage=henonp(Himage,p,q,r,s,e,mysiz);
 p=.205; q=.265; r=.05; s=200;
 Himage=henonp(Himage,p,q,r,s,e,mysiz);
 p=.1; q=.24; r=.05; s=80;
 Himage=henonp(Himage,p,q,r,s,e,mysiz);
 p=0; q=.33; r=.03; s=60;
 Himage=henonp(Himage,p,q,r,s,e,mysiz);
 p=.3; q=0; r=.01; s=60;
 Himage=henonp(Himage,p,q,r,s,e,mysiz);
 image(Himage);
 text(50,50,strcat('energy=',num2str(e)));
 input('next screen','s');
 hold off
 clf
 hold on
 Himage = ones(mysiz,mysiz,3);
 e = 1/8;
 p=.66; q=0; r=.05; s=100;
 Himage=henonp(Himage,p,q,r,s,e,mysiz);
 p=.5; q=0; r=.03; s=200;
 Himage=henonp(Himage,p,q,r,s,e,mysiz);
 p=.205; q=.265; r=.05; s=200;
 Himage=henonp(Himage,p,q,r,s,e,mysiz);
 p=.1; q=.24; r=.05; s=80;
 Himage=henonp(Himage,p,q,r,s,e,mysiz);
 p=0; q=.33; r=.03; s=60;
 Himage=henonp(Himage,p,q,r,s,e,mysiz);
 p=.33; q=0; r=.01; s=60;
 Himage=henonp(Himage,p,q,r,s,e,mysiz);
 image(Himage);
 text(50,50,strcat('energy=',num2str(e)));
 input('next screen','s');
 hold off
 clf
 hold on
 Himage = ones(mysiz,mysiz,3);
 e = 1/6; 
 p=.99; q=0; r=.05; s=100;
 Himage=henonp(Himage,p,q,r,s,e,mysiz);
 p=.66; q=0; r=.04; s=100;
 Himage=henonp(Himage,p,q,r,s,e,mysiz);
 p=.5; q=0; r=.05; s=100;
 Himage=henonp(Himage,p,q,r,s,e,mysiz);
 p=.205; q=.265; r=.05; s=200;
 Himage=henonp(Himage,p,q,r,s,e,mysiz);
 p=.1; q=.24; r=.05; s=80;
 Himage=henonp(Himage,p,q,r,s,e,mysiz);
 p=0; q=.33; r=.02; s=60;
 Himage=henonp(Himage,p,q,r,s,e,mysiz);
 p=.3; q=0; r=.01; s=60;
 Himage=henonp(Himage,p,q,r,s,e,mysiz);
 image(Himage);
 text(50,50,strcat('energy=',num2str(e)));
 hold off
 
function [Himage]=henonp(Himage,y,y1,d,iter,e,mysiz)
range=1;
x=0;
x1=real(sqrt(2*e-y1*y1-y*y+2*y*y*y/3));
i=1;
while i<=iter*100
    xh=x;
    x2=-x-2*x*y;
    y2=-y-x*x+y*y;
    x1=x1+x2*d;
    y1=y1+y2*d;
    x=x+x1*d;
    y=y+y1*d;
    if x*xh<=0
        i=i+1;
        if (abs(y)<range)&&(abs(y1)<range)
            Himage(int16(floor(mysiz*(y+range)/2/range)),int16(floor(mysiz*(y1+range)/2/range)),:)=0;    
        end
    end
end