function  dflDE(A,a,m,e,title)

% This function m-file plots the direction field of a 2 by 2 by linear system 
% of DEs of the form  x'(t)=A x(t), where x in R^2 and A is real 2 by 2.
% 
% Inputs: 2 by 2 Matrix A; real number a which gives the boundary of the plot:
%         -a <= x, y <= a  in uniform steps of size  h= (2a)/(n-1); 
%          m = # nodes per line (square grid) (suggestion: set m = 9, 13 or 21)
%         If e = 1: we draw the eigenspaces for real eigenvalues and for complex 
%         eigenvalues, we draw the major axes of the elliptical spirals from an
%         orthogonal set of real and imaginary eigenvector parts into the plot.
%         (e is optional)
%         If title = 0, we plot without a title line. The default is set to 
%         title = 1.  (Specifying 'title' is optional) 
% Output: Direction field plot; input data for A.        

if nargin < 4, e = 0; end, if nargin < 5, title = 1; end,
[n n] = size(A);                                              % initialize

if n ~= 2,  'A is not 2 by 2',  return, end;                  % data checks
if a <= 0
 'right interval endpoint is negative; bad choice; a > 0 only', return, end  
hold off                                                      % clear screen

h = 2*a/(m-1); edge=(a+h)/a; x = -a:h:a; y = x; len = 0.6*h;  % initialize grid
plot([0,0],edge*[-a,a],'-'); axis(edge*[-a a -a a]), hold on; % hold images
axis('square'); axis auto; plot(edge*[-a,a],[0,0],'-');       % axis plotting
xlabel('x2'); ylabel('x1','Rotation', 0);                     % labels
if title == 1, if a == 5                                      % for a = 5 :
text(-1.32*a, 1.35*a,' Direction Field for the Linear System of DEs '), 
text(0.05*a, 1.35*a,' d(x(t))/dx = A * x(t);  A =  ') % adjust name of matrix
text(0.85*a,1.39*a,[num2str(A(1,1)),'   ', num2str(A(1,2))]);
text(0.85*a,1.31*a,[num2str(A(2,1)),'   ', num2str(A(2,2))]); end
if a ~= 5                                                     % for a ~= 5 :
text(-1.35*a, 1.39*a,' Direction Field for the Linear System of DEs ') 
text(0.06*a, 1.39*a,' d(x(t))/dx = A * x(t);  A =  ') % adjust name of matrix
text(0.88*a,1.43*a,[num2str(A(1,1)),'   ', num2str(A(1,2))]);
text(0.88*a,1.35*a,[num2str(A(2,1)),'   ', num2str(A(2,2))]); end, end

for i=1:m, plot(x(i)*ones(1,m),y,'+'); end                    % plot mesh points

for i=1:m                              % plot scaled tangent vectors
  for j=1:m      
    t = A*[x(i);y(j)];                                % tangent vector x'(t)
    if norm(t) ~= 0, t = len*t/norm(t); end           % normalized
    if norm(t) ~= 0
       plot([x(i);x(i)+t(1)],[y(j);y(j)+t(2)],'-');  end, end, end
  
% follow with eigenvector plots, if desired:

if e == 1
  [v d]=eig(A);  
  if norm(d-d') == 0                                  % for real eigenvalues
    v1 = v(:,1); E1=1.1*a*v1/norm(v1,inf); 
    plot([E1(1),-E1(1)],[E1(2),-E1(2)],':');
    v2 = v(:,2); E2=1.1*a*v2/norm(v2,inf); 
    plot([E2(1),-E2(1)],[E2(2),-E2(2)],':'); 
  else                                                % for complex eigenvalues
    vr = real(v(:,1)); vi = imag(v(:,1)); vrdotvi = vr'*vi;
    if vrdotvi ~= 0                                   % find principal axes
      whirl = -(vr'*vr-vi'*vi)/(2*vrdotvi);
      whirl = whirl + sqrt(whirl^2 + 1);
      VR = whirl*vr - vi; vi = whirl*vi + vr; vr = VR; end
    E1=1.1*a*vr/norm(vr,inf);   
    plot([E1(1),-E1(1)],[E1(2),-E1(2)],':');
    text(1.05*E1(1),1.05*E1(2),'u');
    E2=1.1*a*vi/norm(vi,inf);
    plot([E2(1),-E2(1)],[E2(2),-E2(2)],':'); 
    text(1.05*E2(1),1.05*E2(2),'w');
    end , end