[x, y] = meshgrid(-3:.2:1, -1.5:.2:3);
z = y.^2 - x.^2 + x;

contour(x, y, z, 10)

[u, v] = gradient(z, 0.1);
hold on
quiver(x,y, u, v)

xspan = 0:.1:5;
y0 = 3;

[xx,yy] = ode45(@ODE_eq, xspan, y0);

plot(xx, yy, 'm','LineWidth',3)


function Dy = ODE_eq(x,y)
    Dy(1) = -sin(y.*x);
end