%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Math 3520: Differential Equations and Dynamical Systems
%   Nonlinear first-order systems: X'=F(X)
%   Flow fields for Example 6.1.1: x'=x+exp(-y), y'=-y
%   Trajectories using Runge-Kuta 4th order solver
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



%%%% Runge-Kutta 4th order solver:
function X = rk4a(X0,F,steps,h)
 % Input
 %  X0 = [x0,y0]   Row vector, initial position
 %  F  = function handle, computes Xdot(x,y)
 %  steps = number of time steps to compute (e.g., 100)
 %  h = time increment per step (e.g., 1/steps )
 % Output
 %  X = (1+steps)x(2) matrix of (x,y) values along the
 %                    trajectory starting at X0
  X = zeros(1+steps,2);  % Solution to be computed
  X(1,:) = X0;  % Initial point
  for n=1:steps
    f1 = F(X(n,:));
    f2 = F(X(n,:)+f1*h/2);
    f3 = F(X(n,:)+f2*h/2);
    f4 = F(X(n,:)+f3*h);
    X(n+1,:) = X(n,:) +h*(f1+2*f2+2*f3+f4)/6;
  end
end

%%%% Example 6.1.1 from Strogatz, 3rd edition, p.163
Xdot = @(V) [V(1)+exp(-V(2)), -V(2)];

%%%% Flow field 
% Initialize the grid in phase space where arrows will
% indicate the flow X'=(x',y') at position X=(x,y):
[x,y]=meshgrid(-4:4,-4:4);  % 2-D grid, integers, 9x9 
dx = x+exp(-y); dy = -y;  % Raw differentials
dd = sqrt(dx.*dx + dy.*dy) + 0.01;  % Raw speeds, nonzero 
dy = dy./dd; dx = dx./dd; % Normalized unit arrows
ral = 0.4;                 % Relative arrow length
quiver(x,y,dx,dy,ral);

%%%% Nullclines
% Choose x,y ranges to plot nullclines within the flow field
xnc=-4:0.1:4; ync=-log(4):0.1:4;  % to stay in [-4,4]x[-4,4]
hold on    % to superimpose on the flow field graph
plot(xnc,0*xnc,"--")        % y=0 nullcline, for -4<x<4
plot(-exp(-ync),ync,"--")   % x=-exp(-y) nullcline, -log(4)<y<4

%%%% Trajectories
hold on;  % Plot all trajectories atop the flow field
h=0.01;   % Use same time step for same accuracy
X = rk4a([-3,4], Xdot,  50,h); plot(X(:,1),X(:,2))
X = rk4a([-2,4], Xdot, 100,h); plot(X(:,1),X(:,2))
X = rk4a([-1,4], Xdot, 200,h); plot(X(:,1),X(:,2))
X = rk4a([-0.5,4], Xdot, 250,h); plot(X(:,1),X(:,2))
X = rk4a([ 0.5,4], Xdot, 200,h); plot(X(:,1),X(:,2))
X = rk4a([ 0,4], Xdot, 300,h); plot(X(:,1),X(:,2))
X = rk4a([ 1,4], Xdot, 150,h); plot(X(:,1),X(:,2))
X = rk4a([ 2,4], Xdot,  90,h); plot(X(:,1),X(:,2))
X = rk4a([ 3,4], Xdot,  60,h); plot(X(:,1),X(:,2))

X = rk4a([-1.1,0], Xdot, 380,h); plot(X(:,1),X(:,2))
X = rk4a([-0.9,0], Xdot, 410,h); plot(X(:,1),X(:,2))

X = rk4a([-4,-2], Xdot, 150,h); plot(X(:,1),X(:,2))
X = rk4a([-4,-2.2], Xdot, 250,h); plot(X(:,1),X(:,2))
X = rk4a([-4,-2.5], Xdot, 250,h); plot(X(:,1),X(:,2))
X = rk4a([-4,-3], Xdot, 100,h); plot(X(:,1),X(:,2))
X = rk4a([-4,-4], Xdot,  25,h); plot(X(:,1),X(:,2))


title("Flow Field and Trajectories for Example 6.1.1")