% pplane2d.m 
% Plots a phase plane for any 2D ODE system that has been defined via 
% INITSYS.M 

cla 
F=eqn_name; 
ng=26; 
ax=initax; 
hold on 
x1g=linspace(ax(1),ax(2),ng); 
x2g=linspace(ax(3),ax(4),ng); 
dx1=x1g(2)-x1g(1); 
dx2=x2g(2)-x2g(1); 
x=[x1g(1); x2g(1)]; 
t=0; 

for j=1:ng 
 x(1)=x1g(j); 
 for k=1:ng 
  x(2)=x2g(k); 
  V=feval(F,t,x); 
  v=V(2)/V(1); 
  d=.3/max(1/dx1,abs(v)/dx2); 
  plot([x(1)-d,x(1)+d],[x(2)-v*d,x(2)+v*d],'g') 
 end 
end 
axis(ax) 
title(['Phase Plane, ODE45,  ',F,'.m']) 
xlabel('x2 (vertical) vs x1 (horizontal)')