function phaseportrait;

% This function plots phase portraits.  Users need to modify the
% following:
% (1)  Bounds for the plot need to be "intelligently" specified
%       to include all relevant features that one wants to display
%       This usually means that some analytical information needs
%       to be known in advance.
% (2)  The vector field of interest needs to be specified in the
%        subfunction called "xdot".  It is essential that one uses
%        MATLAB array operations, i.e. .*  ./  .^ instead of * / ^
% Of course, many customizations are possible.

xmin = -2;   %bounds for plotting
xmax = 2;
ymin = -2;
ymax = 2;

x = [xmin:(xmax-xmin)/30:xmax];
y = [ymin:(ymax-ymin)/30:ymax];

[X,Y] = meshgrid(x,y);  %creates matricies full of x,y coordinates
[m,n] = size(X);    %find dimensions of matrix X

%define the vector field using the function f(t,x)
Z = f(0 , [reshape(X,1,m*n);  reshape(Y,1,m*n)]);
U = Z(1,:);     %i component
U = reshape(U,m,n);
V = Z(2,:);       %j component
V = reshape(V,m,n);

%make into direction field by rescaling
U = U./sqrt(U.^2 + V.^2 + 1.0e-14);
V = V./sqrt(U.^2 + V.^2 + 1.0e-14);
%plot  vector field
quiver(X,Y,U,V);
axis([xmin xmax ymin ymax]);
hold on;

%plot sample trajectories

options = odeset('events',@checkbounds);
for x0 =xmin:(xmax-xmin)/8:xmax,
   for y0 = ymin:(ymax-ymin)/8:ymax,
       [t,z,d1,d2,d3] = ode45(@f,[0 100],[x0 y0],options);   % solve with (x0,y0) as initial data
       plot(z(:,1),z(:,2),'k');    %plot solution
   end;
end;
hold off;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%define the vector field (remember to use .'s for array operations!)

function xdot = f(t,x);

X = x(1,:);  %first row has 'x' values
Y = x(2,:);  %second row has 'y' values

U = 2*X.*Y;           %here's where the vector field is defined
V =  Y.^2 - X.^2; 

xdot = [U; V];    %first column is i component, second j component

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  this function allows ode45 to terminate when (x,y) gets too far out of bounds
function [value,isterminal,direction] = checkbounds(t,x);
X = x(1,:);  %first row has 'x' values
Y = x(2,:);  %second row has 'y' values
value(1) = 20 - sqrt(X^2+Y^2);   %if more than certain distance from origin, terminate
isterminal(1) = 1;   %terminate if value == 0
direction(1) = 0;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%