<< Graphics`PlotField`

<< Calculus`VectorAnalysis`

SetCoordinates[Cartesian[x, y, z]];

track[s_] := {s, -Sin[s]}

potential[x_, 
        y_] := -(9.8*y + a (((Sin[x] + y))^2))

FullSimplify[Grad[potential[x, y]]]

accel[x_, 
        y_] := {-2 a Cos[x] (y + Sin[x]), -9.8` - 
          2 a (y + Sin[x]), 0}

accx[x_, y_] := accel[x, y] . {1, 0, 0}

accy[x_, y_] := accel[x, y] . {0, 1, 0}

g1 = 
      PlotGradientField[
        potential[x, y], {x, -2 π, 2  π}, {y, -4, 
          4}]

g2 = 
      ParametricPlot[track[s], {s, -2 π, 2  π}, 
        Compiled -> False]

Show[g1, g2]

a = 50;

k = 4.6; dk =  .05;

solution = 
      NDSolve[{x''[t] == accx[x[t], y[t]], y''[t] == 
            accy[x[t], y[t]], x'[0] == 0, y'[0] == 0, 
          y[0] == 0, x[0] == 0}, {x, y}, {t, 0, k}]

g4 = 
      ParametricPlot[{x[t] /. solution[[1]], 
          y[t] /. solution[[1]]}, {t, 0, k}, 
        Compiled -> False]

MatrixForm[
      Table[{x[a] /. solution[[1]], y[a] /. solution[[1]]}, {a, dk, 
          k, dk}]]

A = Join[{{0, 0}}, %];

range = {{Min[Table[A[[a, 1]], {a, 1, k/dk + 1}]] -  .1, 
          Max[Table[A[[a, 1]], {a, 1, k/dk + 1}]] + 
            1}, {Min[Table[A[[a, 2]], {a, 1, k/dk + 1}]] -  .1, 
          Max[Table[A[[a, 2]], {a, 1, k/dk + 1}]] + 1}}

g1 = 
      ParametricPlot[{track[s]}, {s, 0, π}, AspectRatio -> Automatic, 
        Compiled -> False]

ext = 
      Table[Show[{g1, Graphics[{PointSize[ .05], Point[A[[t]]]}]}, 
          AspectRatio -> Automatic, PlotRange -> range], {t, 1, 
          k/dk}]

Export["Loose_Parabola5.GIF", %, 
      AnimationDisplay -> dk]