Clear[a, b, x, y, t, endtime];

surfaceplot = Plot3D[f[x, y], {x, -2, 2}, {y, -2, 2}, PlotRange -> All, BoxRatios -> Automatic, MaxRecursion -> 4, AxesLabel -> {"x", "y", "z"}];

{a, b} = {.5, -1};
startingpoint3Dplot = Graphics3D[{Red, Sphere[{a, b, f[a, b]}, 0.2]}];
gradf[x_, y_] = {D[f[x, y], x], D[f[x, y], y]};
equationx = x'[t] == gradf[x[t], y[t]][[1]];
equationy = y'[t] == gradf[x[t], y[t]][[2]];
starterx = x[0] == a;
startery = y[0] == b;
endtime = 80;
approxsolutions = NDSolve[{equationx, equationy, starterx, startery}, {x[t], y[t]}, {t, 0, endtime}];
gradtrajectory[t_] = {x[t] /. approxsolutions[[1]], y[t] /. approxsolutions[[1]]};
gradtrajectory3D[t_] = Join[gradtrajectory[t], {f @@ gradtrajectory[t]}];
gradtrajectory3Dplot = ParametricPlot3D[gradtrajectory3D[t], {t, 0, endtime}, PlotStyle -> {Orange, Thickness[0.015]}];
Show[surfaceplot, gradtrajectory3Dplot, startingpoint3Dplot,  PlotRange -> All, AxesLabel -> {"x", "y", "z"}]