clear all
% we open a new figure and set the axis limits
figure(1)
clf
hold on
axis equal

% global phase portrait
axis([-1,2,-1,1])
% local phase portrait around fixed point (0,0)
%axis([-.1,.1,-.1,.1])
% local phase portrait around fixed point (1,0)
%axis([1-.1,1+.1,-.1,.1])
% local phase portrait around fixed point (1/2,1/4)
%axis([.5-.1,.5+.1,.25-.1,.25+.1])

% we define the function f defining the ODE system x'=f(t,x)
% the inaxis function is added to avoid solving the ODE outside
% the current axis limits
f = inaxis(@(t,x) [x(1)^2-x(1)^3-x(1)*x(2);x(1)*x(2)-x(2)/2]);

% we start a loop to compute and draw 100 solutions of the ODE
% with various initial points.
for k=1:100
    % first we pick a point at random 
    x0 = rand(2,1); 
    % this gives a random point in (0,1)^2. 
    % Now we scale it and translate it to get a random point
    % in the window corresponding to the axis limits
    ax = axis;
    x0(1) = x0(1)*(ax(2)-ax(1))+ax(1);
    x0(2) = x0(2)*(ax(4)-ax(3))+ax(3);
    % next we solve the ODE from t=0 to t=15 with this initial point,
    % then plot the solution and plot arrows to show directions
    [T,x] = ode45(f,[0,15],x0);
    plot(x(:,1),x(:,2))
    plotarrows(x)
    % we repeat the operation solving the ODE backward from t=0 to t=-15
    [T,x] = ode45(f,[0,-15],x0);
    plot(x(:,1),x(:,2))
    % here we have to flip the sequence of points, otherwise the arrows
    % will correspond to the backward direction.
    x = flipud(x);
    plotarrows(x)
end

% finally we add to the plot the nullclines for this example
x1 = linspace(-1,2,100);
plot(x1,x1-x1.^2,'r','Linewidth',3)
plot([0,0],[-1,1],'r','Linewidth',3)
plot([.5,.5],[-1,1],'g','Linewidth',3)
plot([-1,2],[0,0],'g','Linewidth',3)



