In[]:=
k=266*1000;m=3755;r=0.5;L=2;A=L^2/4+L*x[t]+x[t]^2+y[t]^2-0.5*L^2*Cos[theta[t]]-L*x[t]*Cos[theta[t]]+0.25*L^2*Cos[theta[t]]^2-L*y[t]*Sin[theta[t]]+0.25*L^2*Sin[theta[t]]^2;B=L^2/4-L*x[t]+x[t]^2+y[t]^2-0.5*L^2*Cos[theta[t]]+L*x[t]*Cos[theta[t]]+0.25*L^2*Cos[theta[t]]^2+L*y[t]*Sin[theta[t]]+0.25*L^2*Sin[theta[t]]^2;oneA=0.5*k*L^2*Sin[theta[t]]+0.25*k*L^2*(-2*Cos[theta[t]]*Sin[theta[t]])+0.25*k*L^2*(2*Sin[theta[t]]*Cos[theta[t]]);oneB=-k*r*(0.5*L^2*Sin[theta[t]]+L*x[t]*Sin[theta[t]]+0.25*L^2*(-2*Cos[theta[t]]*Sin[theta[t]])-L*y[t]*Cos[theta[t]]+0.25*L^2*(2*Sin[theta[t]]*Cos[theta[t]]))/(2*Sqrt[A]);oneC=-k*r*(0.5*L^2*Sin[theta[t]]-L*x[t]*Sin[theta[t]]+0.25*L^2*(2*Cos[theta[t]]*Sin[theta[t]])+L*y[t]*Cos[theta[t]]+0.25*L^2*(2*Sin[theta[t]]*Cos[theta[t]]))/(2*Sqrt[B]);ONE=oneA+oneB+oneC;two=2*k*x[t]-k*r*(L+2*x[t]-L*Cos[theta[t]])/(2*Sqrt[A])-k*r*(-L+2*x[t]+L*Cos[theta[t]])/(2*Sqrt[B]);three=-9.81*m+2*k*y[t]-k*r*(2*y[t]-L*Sin[theta[t]])/(2*Sqrt[A])-k*r*(2*y[t]+L*Sin[theta[t]])/(2*Sqrt[B]);S1=NDSolve[{x''[t]==-two/(m),y''[t]==-three/(m),theta''[t]==-ONE*4/(m*L^2),x[0]==0.0522,y[0]==0.6,theta[0]==0.0349,x'[0]==0,y'[0]==0,theta'[0]==0},{x,y,theta},{t,0,30},Method->{"TimeIntegration"->{"ExplicitRungeKutta","DifferenceOrder"->8}}];S2=NDSolve[{x''[t]==-two/(m),y''[t]==-three/(m),theta''[t]==-ONE*4/(m*L^2),x[0]==0.0522,y[0]==0.6,theta[0]==0.0349,x'[0]==0,y'[0]==0,theta'[0]==0},{x,y,theta},{t,0,30},Method->{"TimeIntegration"->"ExplicitEuler"}];S3=NDSolve[{x''[t]==-two/(m),y''[t]==-three/(m),theta''[t]==-ONE*4/(m*L^2),x[0]==0.0522,y[0]==0.6,theta[0]==0.0349,x'[0]==0,y'[0]==0,theta'[0]==0},{x,y,theta},{t,0,30}];X1=Evaluate[{x[t]}/.S1];X2=Evaluate[{x[t]}/.S2];X3=Evaluate[{x[t]}/.S3];Plot[{X1,X2,X3},{t,0,30},PlotLegends->{"Explicit Runge Kutta","Explicit Euler","Mathematica optimised"}]
Out[]=
