f[x_, y_] = E^((1/16) (x^2 + 4 y^2))
m[x_, y_] = D[f[x, y], x];
n[x_, y_] = D[f[x, y], y];
{a, b} = {-.5, -.5};
starterpoint = {a, b};
Clear[x, y, t];
equationx = x'[t] == m[x[t], y[t]];
equationy = y'[t] == n[x[t], y[t]];
starterx = x[0] == a;
startery = y[0] == b;
endtime = 2;
approxsolutions = NDSolve[{equationx, equationy, starterx, startery}, {x[t], y[t]}, {t, 0, endtime}];
Clear[trajectory];
trajectory[t_] = {x[t] /. approxsolutions[[1]], y[t] /. approxsolutions[[1]]};
trajectoryplot = ParametricPlot[trajectory[t], {t, 0, endtime}, PlotStyle -> {{Red, Thickness[0.015]}}];
starterplot = Graphics[{Red, PointSize[0.06], Point[starterpoint]}];
Show[gradfieldplot, starterplot, trajectoryplot, PlotRange -> All, Axes -> True, AxesLabel -> {"x", "y"}]


{x[t_], y[t_]} = {a, b} + t gradf @@ {a, b}
surfaceplot = Plot3D[f[x, y], {x, -1, 1}, {y, -1, 1}];
traj = ParametricPlot3D[{x[t], y[t], f[x[t], y[t]]}, {t, 0, 2}, PlotStyle -> {Red, Thickness[0.01]}];
Show[traj, surfaceplot, PlotRange -> Automatic]