How to add direction line into the phase plane plot?

Yes, how to add direction line into the phase plane plot? I have the codes below which it work successful but I want to add some direction arrow into my phase plane.

function phase()
IC = [1 1;1 2;1 3;1 4;1 5];
hold on
for ii = 1:length(IC(:,1))
    [~,X] = ode45(@EOM,[-5 5],IC(ii,:));
    u = X(:,1);
    w = X(:,2);
    plot(u,w,'r')
end
xlabel('u')
ylabel('w')
grid
end
function dX = EOM(t, y)
dX = zeros(2,1);
u  = y(1);
w  = y(2);
A  = 1;
B  = 1;
dX = [(w*u^2 - B*u);...
      (A - w - w*u^2)];
end

And the plot below is from the codes above:

And I want to add the direction line, which something like this:

%

Answer 1

[Mischa Kim]

Alex, use quiver. In

 quiver(x,y,u,v)

the x, y are the grid points, u and v you get from the dX vector from the differential equations function.

[Alex]

Where should I put quiver(x,y,u,v) into? At the begining of the function or before plot?

[Alex]

 Something like tihs?

function dX = phase(t, y);
  u  = y(1);
  w  = y(2);
  A  = 1;
  B  = 1;
  dX = [(w*u^2 - B*u);...
        (A - w - w*u^2)];
  end
  dX = @(t,Y) [Y(2); -sin(Y(1))];
  y1 = linspace(-2,8,20);
  y2 = linspace(-2,2,20);
  [x,y] = meshgrid(y1,y2);
  size(x)
  size(y)
  u1 = zeros(size(x));
  w1 = zeros(size(x));
  % we can use a single loop over each element to compute the derivatives at
  % each point (y1, y2)
  t=0; % we want the derivatives at each point at t=0, i.e. the starting time
  for i = 1:numel(x)
      Yprime = f(t,[x(i); y(i)]);
      u1(i) = Yprime(1);
      w1(i) = Yprime(2);
  end
  quiver(x,y,u,w,'r'); figure(gcf)
  xlabel('y_1')
  ylabel('y_2')
  axis tight equal;

It doesnt work properly tho!

[Mischa Kim]

Using the same EOM function I'd do something like

 function phase()

 IC = 5*(rand(50,2)-0.5);
 hold on
 for ii = 1:length(IC(:,1))
    [~,X] = ode45(@EOM,[-5 5],IC(ii,:));
    u = X(:,1);
    w = X(:,2);
    plot(u,w,'r')
 end
 xlabel('u')
 ylabel('w')
 grid

 x       = -4:0.5:4;
 y       = -4:0.5:4;
 [xg,yg] = meshgrid(x,y);
 dxg     = yg.*xg.^2 - xg;
 dyg     = ones(length(x)) - yg - yg.*xg.^2;
 scale   = 5;
 quiver(xg,yg,dxg,dyg,scale)

 end

[Alex]

Thanks again, you are the best !!