function PhasePortrait(dydt,y1,y2,s,Fmx)

% 2D phase portrait for
%
%   [dy(1)/dt,dy(2)/dt] = dydt(y(1),y(2))
% 
% Inputs; 
%  dydt     derivative function
%  y1, y2   vectors of values defining all coordinate pairs (y1,y2)
%           at which to plot the vector F, scaled based on s.
%  s        Scaling ratio of max vector length to mesh spacing (try 0.5)
%  Fmx      Vectors F with |F|>Fmx will get scaled to magnitude Fmx
%           If Fmx=0, no thresholding on |F| is performed.
% 
% Orbits, axis labels, and title can be added by later plotting calls.  

  [Y1,Y2]=meshgrid(y1,y2);
  siz = size(Y1);
  t = 0;  
  F = feval(dydt,t,[(Y1(:))'; (Y2(:))']);
  if (Fmx>0) 
    Fnorm = sqrt(F(1,:).^2 + F(2,:).^2);
    Fnorm(Fnorm<1e-6*Fmx) = 1e-6*Fmx;
    notbig = Fnorm<Fmx;
    Fnorm(notbig) = Fmx;
    F1 = reshape(F(1,:)./Fnorm,siz);
    F2 = reshape(F(2,:)./Fnorm,siz);
    quiver(Y1(~notbig),Y2(~notbig),F1(~notbig),F2(~notbig),s,'m')
    hold on
    quiver(Y1(notbig),Y2(notbig),F1(notbig),F2(notbig),s,'b')
    hold off
  else
    F1 = reshape(F(1,:),siz);
    F2 = reshape(F(2,:),siz);
    quiver(Y1,Y2,F1,F2,s)
  end