%
%  Linear 2D ODE system example
%  Demonstrates use of phase space and vector fields
%  Solve x' = -y, y' = -x, with initial condition x(0)=2, y(0)=1.8
%  Victor Matveev, October 9, 2015
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear;                                % Erase everything
figure('position',[500 500 800 300]); % Optional: specify window size

f = @(x,y) -y;                        % Define x'=f(x,y) and ...
g = @(x,y) -x;                        %   ...  y'=g(x,y) velocity functions

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%              Plot vector field (flow field) using "quiver"
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

subplot(2,2,[1 3]); hold on;         % Next plot will be in 2x2 "cells" #1 and #3
u = -2 : 0.5 : 2;                    % Grid of points to use along x- and y-axes
[X, Y] = meshgrid(u, u);             % Create X & Y coordinate mesh covering a grid
quiver(X, Y, f(X, Y), g(X, Y));      % <<== Plot the vector field: easy!
xlabel('x'); ylabel('y');            % Add nice labels and title
title('Vector field F(x,y) = (-y, -x)');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%       Use Euler Method to find numerical (approximate) solution
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

tTotal = 3; h = 0.02; N = round(tTotal / h); % Usual Euler parameters
h = tTotal / N;                              % Re-calculate h, just in case
x(1) = 2;  y(1) = 1.8;  t(1) = 0;            % The initial conditions (see above)

for n = 2 : N + 1
    x(n) = x(n-1) + h * f(x(n-1), y(n-1));   % No more complicated than...
    y(n) = y(n-1) + h * g(x(n-1), y(n-1));   %   ...solving a 1D problem!
    t(n) = h*(n - 1);                        % Can move this outside if we like...
    
    subplot(2,2,[1 3]);               % Plot the trajectory in cells #1 and #3
    plot(x, y, 'r-', 'linewidth', 2); % Plot trajectory (x(t), y(t)) on phase plane
    axis tight;                       % Make coordinate ranges as narrow as possible

    subplot(2,2,2);                   % Next plot will go to cell #2 of 2x2 grid
    plot(t, x, 'r-', 'linewidth', 2); % Plot x(t) using red line of width 2
    xlabel('time'); ylabel('x(t)');   % Add nice labels and title
    title('Solution x(t) as a function of time');
    axis([0 tTotal 0 2]);             % Specify coordinate ranges
    
    subplot(2,2,4);                   % Next plot will go to cell #4 of 2x2 grid
    plot(t, y, 'r-', 'linewidth', 2); % Plot y(t) using red line of width 2
    xlabel('time'); ylabel('y(t)');   % Add nice labels and title
    title('Solution y(t) as a function of time');
    axis([0 tTotal -2 2]);            % Specify coordinate ranges
    
    drawnow;
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%    Calculate and plot the exact solution (with circles), using
%   Linear algebra result: Y(t) = C1 V1 exp(L1*t) + C2 V2 exp(L2*t)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear Y;                     % Erase Y so that we can re-use it again
T = linspace(0, tTotal, 20); % Create 20 data points between 0 and tTotal
V1 = [1 1];   L1 = -1;       % First  eigenvector-eigenvalue pair
V2 = [1 -1];  L2 = 1;        % Second eigenvector-eigenvalue pair
C1 = 1.9;     C2 = 0.1;      % ...So that C1*V1 + C2*V2 = (2, 1.8) = initial condition

for n = 1 : length(T)        % We could vectorize this, if we wanted to...
    Y = C1 * V1 * exp(L1*T(n)) + C2 * V2 * exp(L2*T(n));
    
    subplot(2,2,[1 3]); hold on;
    plot(Y(1), Y(2), 'ok');  % Place a circle on the trajectory
    subplot(2,2,2); hold on;
    plot(T(n), Y(1), 'ok');  % Place a circle on the x(t) plot
    subplot(2,2,4); hold on;
    plot(T(n), Y(2), 'ok');  % Place a circle on the y(t) plot
   
    drawnow;
end;

legend('Euler', 'Exact');  % Add a legend on one of the plots
