step16;

h=0.0001;
fdu=(f(1+h*u(1),1+h*u(2))-f(1,1))/h;
quiver3(1,1,f(1,1),u(1),u(2),fdu,2,'color','k','linewidth',3,'maxheadsize',2);
