function steep_ascent(x,y,z,px,py,N,xstart,ystart)

xsa(1)= xstart;
ysa(1)= ystart;
zsa(1)=func(xsa(1),ysa(1));

for i=1:N
  gradx = ( func(xsa(i)+eps,ysa(i))-func(xsa(i),ysa(i)) )/eps;
  grady = ( func(xsa(i),ysa(i)+eps)-func(xsa(i),ysa(i)) )/eps;
 

  delta = fmin('fgrad',-10,10,[0,1.e-4],xsa(i),ysa(i),gradx,grady);
  xsa(i+1) = xsa(i) + delta*gradx;
  ysa(i+1) = ysa(i) + delta*grady;
  zsa(i+1) = func(xsa(i+1),ysa(i+1));
end

hold off
contour(x,y,z,10)
hold on
quiver(x,y,px,py)
hold on
plot(xsa,ysa)
S = sprintf('Steepest Ascent: N = %d',N);
title(S)
xlabel('x axis')
ylabel('yaxis')

