function r=phase_portrait(F,X,Y,traj,tau,time)
% Filip Piekniewski 2009
% Plots the vector field, nullclines and sample trajectories.
% F - function handle that computes the vector fiels (see sample fields on www.mat.umk.pl/~philip)
% X, Y - ranges
% traj - number of trajectories that will be computed (originating in random points)
% tau - timestep for computing the trajectories
% time - length of the trajectory (time of integration)
% sample call
% phase_portrait('sample_field4',-4:0.2:4,-4:0.2:4,200,0.02,30)
%
% Work on Matlab and Octave
%
xl=min(X(:));
xu=max(X(:));
yl=min(Y(:));
yu=max(Y(:));

[X,Y]=meshgrid(X,Y);
[U,V]=feval(F,X,Y);

figure;
hold on;
contour(X,Y,U,[0 0],'g','LineWidth',1.5)
contour(X,Y,V,[0 0],'r','LineWidth',1.5)
quiver(X,Y,U,V);
for i=1:traj
  x0=((xu-xl)*rand(1,1))+xl;
  y0=((yu-yl)*rand(1,1))+yl;
  tr=[];
  tr(1,:)=[x0,y0]; 
  for t=0:tau:time
    [dx,dy]=feval(F,tr(end,1),tr(end,2));
    tr(end+1,:)=tr(end,:)+tau*[dx,dy];
    if (tr(end,1)>xu)||(tr(end,1)<xl)||(tr(end,2)>yu)||(tr(end,2)<yl) 
      tr=tr(1:end-1,:);
      break;
    end; 
  end;
  plot(tr(:,1),tr(:,2),'k');
end;
axis([xl xu yl yu]);


