(* Az egyenletben szereplő paramétereket rögzítjük. *)

T=2/15

g=9.81

l=2

(* A fölfüggesztési pontra ható erő hatására létrejövő gyorsulás egy lépcsősfüggvény. Ezt az alábbi módon készítjük el. *)

tred[t_]:=Mod[t,2*T]

A[t_]:=500*(-UnitStep[tred[t]-T]+0.5)

Plot[A[t],{t,0,3}]

(* Megoldatjuk a kezdeti érték problémát. Egyelőre csak a [0,1] intervallumon. *)

sol1=NDSolve[{x'[t]==y[t],y'[t]==(A[t]/l+g/l)*x[t],x[0]==1,y[0]==-0.9},
  {x,y},{t,0,1}];ParametricPlot[Evaluate[{x[t],y[t]}/.sol1],{t,0,1},PlotPoints->100,AxesOrigin->{0,0},AspectRatio->2]

(* Uynazazt a kezdeti érték problémát tekintjük egy 10-szer hosszabb időintervallumon. Reménykeltő: korlátos a trajektória. *)

sol2=NDSolve[{x'[t]==y[t],y'[t]==(A[t]/l+g/l)*x[t],x[0]==1,y[0]==-0.9},
  {x,y},{t,0,10}];ParametricPlot[Evaluate[{x[t],y[t]}/.sol2],{t,0,10},PlotPoints->100,AxesOrigin->{0,0},AspectRatio->2]

(* Egy igen hosszú: [0,100] intervallumon figyelve a trajektóriát azt látjuk, hogy a megoldás korlátos, azaz stabil. *)

sol3=NDSolve[{x'[t]==y[t],y'[t]==(A[t]/l+g/l)*x[t],x[0]==1,y[0]==-0.9},
  {x,y},{t,0,100}];ParametricPlot[Evaluate[{x[t],y[t]}/.sol3],{t,0,100},PlotPoints->100,AxesOrigin->{0,0},AspectRatio->2]