%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  Math 613 * Fall 2018 * Homework 3 * Problem 2 * Victor Matveev
%
%  Note that the system is easily solvable in closed-form, but I use
%  numerical solution for the trajectory, to ensure code re-usability
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure;
L = 1;
x = linspace(-L,L,12);    % Use a sparse set of points for the vector field 
[X, Y] = meshgrid(x, x);  % Create the coordinate grid 
F1 = -X + Y;              % Horizontal velocity 
F2 = -Y;                  % Vertical velocity 

subplot(2,3,[1,2,4,5]); hold on;
quiver(X, Y, F1, F2);      % Plot the vector field
xlabel('x'); ylabel('y');

plot([-L L], [-L, L], 'b-', 'linewidth', 2); % Plot x-Nullcline
plot([-L L], [0   0]  'r-', 'linewidth', 2); % Plot y-Nullcline/eigenvector

f = @(t, y) [-y(1) + y(2); - y(2)];         % Set up the ODE file
[T, Y] = ode45(f, [0 10], [0; 1]);          % Compute the trajectory
plot(Y(:,1), Y(:,2), 'm-', 'linewidth', 2); % Plot the trajectory

legend('Vector field', 'N_x', 'N_y=eigenvector', 'Trajectory');
axis tight;

subplot(2,3,3);
plot(T, Y(:,1), 'm-', 'linewidth', 2);     % Plot x(t)
xlabel('time'); %ylabel('x(t)');
title('x(t)');
axis tight;

subplot(2,3,6);
plot(T, Y(:,2), 'm-', 'linewidth', 2);    % Plot y(t)
xlabel('time'); %ylabel('y(t)');
title('y(t)');
axis tight;