% Modified from (Lynch, 2004):
% Chapter 8 - Planar Systems.
% Program_8b - Phase portrait.
% Copyright Birkhauser 2004. Stephen Lynch.

% Phase portrait of a linear system of ODE's.
% IMPORTANT - Requires vectorfield.m.
clear all
% original --> sys=inline('[x(1)+2*x(2); 3*x(1)+2*x(2)]','t', 'x');  % Define inline function
% sys=inline('[x(2); -x(1)+x(2)*(1-x(1)^2)]','t', 'x');  % van der Pol
 sys=inline('[x(2)+x(2)^2; -x(1)/2 + x(2)/5 - x(1)*x(2) + (6/5)*x(2)^2]','t', 'x');  %  two-eyed monster
options = odeset('RelTol',1e-4,'AbsTol',1e-4);   % set solver tolerances
% original --> [t,xa]=ode45(sys,[0 20],[7.1 .1],options);       % Integrate differential equations
[t,xa]=ode45(sys,[0 100],[7.1 1.1],options);       % Integrate differential equations

fsize=15;                      % set font size for graph
% ----- render x, y vs t graph  -----------
figure(1)
plot(t,xa(:,1),'r')         % plot x vs t, color red
hold on                     % retain the current plot to add more than 1 function
plot(t,xa(:,2),'b')         % plot y vs t, color red
hold off                    % release hold on current plot
legend('x(t)','y(t)')
set(gca,'xtick',[0:5:20],'FontSize',fsize)
set(gca,'ytick',[0:2:8],'FontSize',fsize)
xlabel('time','FontSize',fsize)
ylabel('x, y','FontSize',fsize)

% ----- render x vs y graph with flow field -----------
figure(2)
% warning off MATLAB:divideByZero
vectorfield(sys,-3:.25:3,-3:.25:3)       % draw vector field
    % ----- integrate and plot differential equations -------------
     hold on            
     for x0=-3:1.5:3
         for y0=-3:1.5:3
            [ts,xs] = ode45(sys,[0 5],[x0 y0]);
            plot(xs(:,1),xs(:,2))
         end
     end
     for x0=-3:1.5:3
         for y0=-3:1.5:3
            [ts,xs] = ode45(sys,[0 -5],[x0 y0]);
            plot(xs(:,1),xs(:,2))
         end
     end
     hold off            
    % ------ for a several initial conditions ------------------------
axis([-3 3 -3 3])       % set axis
set(gca,'xtick',[-3:1:3],'FontSize',fsize)
set(gca,'ytick',[-3:1:3],'FontSize',fsize)
xlabel('x(t)','FontSize',fsize)
ylabel('y(t)','FontSize',fsize)
% End of Program.