Henon-Heiles Oscillator

Clear["Global`*"];

eq1 = x''[t] + x[t] + g 2 x[t] y[t] == 0

eq2 = y''[t] + y[t] + g (x[t]^2 - y[t]^2) == 0

g = 1 ; n = 20 ;

eq1

eq2

v = (x^2 + y^2)/2 + g (x^2 y - y^3/3)

Plot[ v /. y->0, {x,-2,2} ] ;

Plot[ v /. x->0, {y,-2,2} ] ;

Plot3D[ v, {x,-2,2},{y,-2,2} ] ;

x0 = 0.3; xp0 = 0; y0 = 0; yp0 = 0.4;

sol = NDSolve[ { eq1, eq2,
                 x[0]==x0, x'[0]==xp0,
                 y[0]==y0, y'[0]==yp0 },
               {x[t],y[t]}, {t,0,n 2Pi},
                MaxSteps->Infinity ] ;

??NDSolve

x1[t] = x[t] /. sol ; y1[t] = y[t] /. sol ;

Plot[ Evaluate[x[t] /. sol], {t,0,n 2Pi} ] ;

Plot[ Evaluate[y[t] /. sol], {t,0,n 2Pi}] ;

Evaluate[ x[t] /. sol] /. t-> N[2Pi]

Evaluate[ y[t] /. sol] /. t-> N[2Pi]

ParametricPlot[ Evaluate[{x[t],y[t]} /. sol],
                {t,0,n 2Pi},
                PlotStyle->RGBColor[0,0,1],
                AspectRatio->Automatic] ;

?ParametricPlot