restart:
f:=(x,y)->0.5*x**2+y:
lit:=plot3d(f(x,y),x=-8..8,y=-50..6,style=patchnogrid):
u:=D(x):v:=D(y):
A:=-(u(t)**2*diff(f(x,y),x,x)(t)+2*u(t)*v(t)*diff(f(x,y),x,y)(t)+v(t)**2*diff(f(x,y),y,y)(t)+1)/(1+diff(f(x,y),x)(t)**2+diff(f(x,y),y)(t)**2):
equa:={D(u)(t)=A*diff(f(x,y),x)(t),
D(v)(t)=A*diff(f(x,y),y)(t),
#position initiale
x(0)=3,y(0)=6,
#vitesse initiale
u(0)=0,v(0)=0}:
solu:=dsolve(equa,{x(t),y(t)},numeric,method=classical,stepsize=0.01,output=listprocedure):
p:=subs(solu,x(t)):
q:=subs(solu,y(t)):
ligne:=spacecurve([seq([p(h),q(h),f(p(h),q(h))],h=seq(0.8*k,k=0..19))],color=blue,thickness=2):
lignecoul1:=pointplot3d([p(0),q(0),f(8,6)],thickness=3,color=red,symbol=circle),seq(pointplot3d([p(h),q(h),f(p(h),q(h))],thickness=3,color=red,symbol=circle),h=seq(0.8*k,k=0..19)):

display([lit,ligne,display(lignecoul1,insequence=true)],style=patchnogrid,scaling=unconstrained,lightmodel=light2,orientation=[-110,85]);