H=0.5*p[t]^2+0.5*q[t]^2;
equations = {q'[t]==D[H,p[t]],p'[t]==-D[H,q[t]]};
init={q[0]==1,p[0]==0};
            

T=2000;
n=T*50;
h=T/n;
sol = NDSolve[{equations,init},{q[t],p[t]},{t,0,T},StartingStepSize->h,MaxSteps->Infinity,
        Method->{"SymplecticPartitionedRungeKutta","DifferenceOrder"->1,"PositionVariables"->{q[t]}}];
		    
ParametricPlot[{q[t],p[t]}/.sol,{t,0,T},PlotRange->{{-1.5,1.5},{-1.5,1.5}}]
