"Implement RK2 for simple pendulum"

(*boundary condition*)
g=9.81;l=1;
theta0=0;omega0=Sqrt[g/l];
theta[0]=theta0;
omega[0]=omega0;
Omega=Sqrt[g/l];
T=2.0*Pi/Omega;
tinit=0;tlast=10T;Nstep=500;
Deltat=(tlast-tinit)/Nstep;
For[i=0,i≤Nstep-1,i++,
thetatudor=theta[i]+(1/2)omega[i]*Deltat;
omegatudor=omega[i]+(1/2)(-g/l)*thetatudor*Deltat;
theta[i+1]=theta[i]+omegatudor*Deltat;
omega[i+1]=omega[i]+(-g/l)*thetatudor*Deltat;
];(*end for i*)

tttheta=Table[{i*Deltat,theta[i]},{i,0,Nstep}];
lptttheta=ListPlot[tttheta,Joined->True,Mesh->All,PlotRange->All,PlotLabel->"theta vs. t"]