% vector and direction field plots

N=[0:0.1:2];
P=[0:0.1:2];

[n,p]=meshgrid(N,P);
% makes a grid of points corresponding to all combinations of N and P
% values
% n and p hold horizontal and vertical components of these points, resp.

dn=n.*(1-p);
dp=p.*(n-1);

% .* does component-wise multiplication, so for each point in the (N,P) 
% grid, we calculate the product N(1-P) and P(N-1)

figure(1)
quiver(n,p,dn,dp)


% normalize velocities
dv=(dn.*dn+dp.*dp).^(1/2);

% again, component-wise operations

dn1=dn./dv;
dp1=dp./dv;
% normalize

figure(2)
quiver(n,p,dn1,dp1)