function Plotslope
x=-5:0.6:5;
y=-5:0.6:5;
delta=0.2;
hold on
axis([-6,6,-6,6])
daspect([1,1,1])
for k=1:length(x)
for m=1:length(y)
eps=delta/(sqrt(1+ff(x(k),y(m))^2));
plot([x(k)-eps, x(k)+eps],...
[y(m)-eps*ff(x(k),y(m)),...
y(m)+eps*ff(x(k),y(m))],'k');
plot(x(k),y(m),'k.','LineWidth',0.2)
end
end
[x0,y0]=ginput(1);
plot(x0,y0,'bo')
[T,Y]=ode45(@ff,[x0,5],y0);
[T1,Y1]=ode45(@ff,[x0,-5],y0);
plot(T,Y,'r',T1,Y1,'r')
function z=ff(x,y)
z=(x+1)*(y+2)*sin(x^2+2*x+pi/3);
end
end