(**************************************)
(*Initialisation*)(* 1D SHM with Euler method *)
Deltat=0.01;tini=0;l=1.0;g=9.81;
omega[0]=0;theta[0]=2*(180/Pi);
T=(2Pi)Sqrt[l/g];
tlast=20T;
(**************************************)

(**************************************)
(* Calculation *)
ksq=g/l;
Ni=IntegerPart[(tlast-tini)/Deltat];
For[i=0,i≤Ni,i++,
t[i]=tini+i*Deltat;
omega[i+1]=omega[i]+Deltat*(-ksq*theta[i]);

theta[i+1]=theta[i]+Deltat*omega[i];

];(*end for *)
tttheta= Table[{t[i],theta[i]},{i,0,Ni}];
ttomega= Table[{t[i],omega[i]},{i,0,Ni}];



(**************************************)

(* Output *)
(**************************************)
ListPlot[tttheta,PlotLabel->"theta vs time "]
ListPlot[ttomega,PlotLabel->"omega vs time "]

(**************************************)

(**************************************)
(*Initialisation*)(* 1D SHM with Euler-Cromer method *)
Deltat=0.01;tini=0;l=1.0;g=9.81;
omega[0]=0;theta[0]=2*(180/Pi);
T=(2Pi)Sqrt[l/g];
tlast=20T;
(**************************************)

(**************************************)
(* Calculation *)
ksq=g/l;
Ni=IntegerPart[(tlast-tini)/Deltat];
For[i=0,i≤Ni,i++,
t[i]=tini+i*Deltat;
omega[i+1]=omega[i]+Deltat*(-ksq*theta[i]);

theta[i+1]=theta[i]+Deltat*omega[i+1];

];(*end for *)
tttheta= Table[{t[i],theta[i]},{i,0,Ni}];
ttomega= Table[{t[i],omega[i]},{i,0,Ni}];



(**************************************)

(* Output *)
(**************************************)
ListPlot[tttheta,PlotLabel->"theta vs time "]
ListPlot[ttomega,PlotLabel->"omega vs time "]

(**************************************)