% Direction field for autonomous 2 by 2 systems m-file (script) 
% 3.2.1998 by HA (modified from Bäckström's 1 eq. case
% Requires the system m-file to be written to work for an initial value
% matrix (with 2 rows) instead of a column vector. This form is perfectly
% good for the ode-solvers
%
clear variables; hold off;
xmax=10; ymax=3; delx=xmax/20; dely=ymax/20;   % *** problem dependent
x=-4:delx:xmax; y=-ymax:dely:ymax;
[X,Y]=meshgrid(x,y);
Z=sys5(0,[X(:),Y(:)]');
DX=Z(1,:);DY=Z(2,:);
DX=reshape(DX,size(X));DY=reshape(DY,size(Y));
  L=sqrt(DX.^2+DY.^2); % Initial lengths
  L=L+1e-10;
  DX1=DX./L; DY1=DY./L; % Unit length for all arrows
%figure(1); quiver(X,Y,DX1,DY1,'.'); grid on;  % The period ('.') suppresses
                                               % arrowheads
figure(1); quiver(X,Y,DX1,DY1,'k'); grid on;
 % axis([0 xmax 0 ymax]); title('Direction field');
hold on
init=ginput(1);[T,Y]=ode45('sys5',[0,10],init);plot(Y(:,1),Y(:,2),'g')
axis([5 7 -1 1])
init=ginput(1);[T,Y]=ode45('sys5',[0,20],init);plot(Y(:,1),Y(:,2),'r')

init=ginput(4);[T,Y]=ode45('sys5',[0,20],init);plot(Y(:,1),Y(:,2),'c')

for i=1:5;
  init=ginput(1);plot(init(1),init(2),'o');[T,Y]=ode45('sys5',[0,20],init);plot(Y(:,1),Y(:,2),'r')
end;
