f[x_,y_]=x^2 y/(x^2+y^2)
cp3D=ContourPlot3D[z==f[x,y],{x,-4.1,4},{y,-4,4},{z,-2,2}]

cp=ContourPlot[f[x,y],{x,-4.1,4},{y,-4,4}]

gradf[x_,y_]=Grad[f[x,y],{x,y}]

vp=VectorPlot[gradf[x,y],{x,-4.1,4},{y,-4,4},VectorPoints->Fine]

Show[cp,vp]

x0=4;y0=2.83;(*2.82*)sol=ParametricPlot[Evaluate[{x[t],y[t]}/.NDSolve[{x'[t]==((2 ((x[t])^3) y[t])/(((((x[t])^2)+((y[t])^2)))^2))-((2 x[t] y[t])/(((x[t])^2)+((y[t])^2))),y'[t]==((2 ((x[t])^2) ((y[t])^2))/(((((x[t])^2)+((y[t])^2)))^2))-(((x[t])^2)/(((x[t])^2)+((y[t])^2))),x[0]==x0,y[0]==y0},{x[t],y[t]},{t,0,30}]],{t,0,20},PlotStyle->Thickness[0.01]]
Show[cp,vp,sol];

FallLinie=NDSolve[{x'[t]==((2 ((x[t])^3) y[t])/(((((x[t])^2)+((y[t])^2)))^2))-((2 x[t] y[t])/(((x[t])^2)+((y[t])^2))),y'[t]==((2 ((x[t])^2) ((y[t])^2))/(((((x[t])^2)+((y[t])^2)))^2))-(((x[t])^2)/(((x[t])^2)+((y[t])^2))),
x[0]==x0,y[0]==y0},{x[t],y[t]},{t,0,30}];
X[t_]=x[t]/.FallLinie[[1,1]];
Y[t_]=y[t]/.FallLinie[[1,2]];
pp3D=ParametricPlot3D[{X[t],Y[t],f[X[t],Y[t]]},{t,0,20}];
Show[cp3D,pp3D]