> restart:with(VectorCalculus):with(plots):with(plottools):
> z:=(x,y)->sin(1/2*x^2-1/4*y^2+3)*cos(2*x+1-exp(y)):
> grad:=VectorCalculus[Gradient](z(x,y),[x,y]);
> plot3d(z(x,y),x=-1.2..1.2,y=-1.2..1.2,axes=normal,numpoints=1000);p3d:=%:
> contourplot(z(x,y),x=-1.2..1.2,y=-1.2..1.2,axes=normal,contours=30,numpoints=3000);cont:=%:
> start:=[-1/4,1/3];ptf[0]:=Vector(start):
> steps:=15;
> for i from 0 to steps do:
  print(ptf[i]):
  pt[i]:=Vector([convert(ptf[i],list)[],z(ptf[i][1],ptf[i][2])]):
  dir[i]:=evalf(Normalize(evalVF(grad,ptf[i])));
  par[i]:=ptf[i]+lambda*dir[i];
  lambd[i]:=fsolve(diff(z(par[i][1],par[i][2]),lambda)=0,lambda=0);
  ptf[i+1]:=eval(par[i],lambda=lambd[i]);
  od:i:='i':
> display(cont,'point(convert(ptf[i],list),color=blue)'$'i'=0..steps,'plot([par[i][1],par[i][2],lambda=0..lambd[i]])'$'i'=0..steps);
> display(p3d,'point(convert(pt[i],list),color=blue,symbol=circle,symbolsize=4)'$'i'=0..steps,'spacecurve([par[i][1],par[i][2],z(par[i][1],par[i][2])],lambda=0