%%%% 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

%Problem 3 (Ex.5.2.[3,...,10], p.91)
%  Plot the flow field in phase space and classify the fixed points of
%   the following systems.  If there are real eigenvectors, plot them
%    on the diagram.
%%% DO THIS FIRST:
% Initialize the grid in phase space where arrows will
% indicate the flow X'=(x',y') at position X=(x,y):
t=-2:0.5:2; [x,y]=meshgrid(t,t);  % 2-D grid 
dt=0.1; % Small increment for short arrows
ral = 0.4;                 % Relative arrow length

% 3(a)  $\dot x=y$, $\dot y=-2x-3y$.
A=[0 1; 2 -3]; A, eig(A)
dx = A(1,1)*x*dt + A(1,2)*y*dt;
dy = A(2,1)*x*dt + A(2,2)*y*dt;
dd = sqrt(dx.*dx + dy.*dy) + 0.01;  % Raw speeds, nonzero 
dy = dy./dd; dx = dx./dd; % Normalized unit arrows
quiver(x,y,dx,dy,ral);
title("Problem 3a");

% 3(b)  $\dot x=5x+10y$, $\dot y=-x-y$.
A=[5 10; -1 -1]; A, eig(A)
dx = A(1,1)*x*dt + A(1,2)*y*dt;
dy = A(2,1)*x*dt + A(2,2)*y*dt;
dd = sqrt(dx.*dx + dy.*dy) + 0.01;  % Raw speeds, nonzero 
dy = dy./dd; dx = dx./dd; % Normalized unit arrows
quiver(x,y,dx,dy,ral);
title("Problem 3b");

% 3(c)  $\dot x=3x-4y$, $\dot y=x-y$.
A=[3 -4; 1 -1]; A, [L,V]=eig(A)
dx = A(1,1)*x*dt + A(1,2)*y*dt;
dy = A(2,1)*x*dt + A(2,2)*y*dt;
dd = sqrt(dx.*dx + dy.*dy) + 0.01;  % Raw speeds, nonzero 
dy = dy./dd; dx = dx./dd; % Normalized unit arrows
quiver(x,y,dx,dy,ral);
title("Problem 3c");

% 3(d)  $\dot x=-3x+2y$, $\dot y=x-2y$.
A=[3 2; 1 -2]; A, eig(A)
dx = A(1,1)*x*dt + A(1,2)*y*dt;
dy = A(2,1)*x*dt + A(2,2)*y*dt;
dd = sqrt(dx.*dx + dy.*dy) + 0.01;  % Raw speeds, nonzero 
dy = dy./dd; dx = dx./dd; % Normalized unit arrows
quiver(x,y,dx,dy,ral);
title("Problem 3d");

% 3(e)  $\dot x=5x+2y$, $\dot y=-17x-5y$.
A=[5 2; -17 -5]; A, eig(A)
dx = A(1,1)*x*dt + A(1,2)*y*dt;
dy = A(2,1)*x*dt + A(2,2)*y*dt;
dd = sqrt(dx.*dx + dy.*dy) + 0.01;  % Raw speeds, nonzero 
dy = dy./dd; dx = dx./dd; % Normalized unit arrows
quiver(x,y,dx,dy,ral);
title("Problem 3e");

% 3(f)  $\dot x=-3x+4y$, $\dot y=-2x+3y$.
A=[-3 4; -2 3]; A, eig(A)
dx = A(1,1)*x*dt + A(1,2)*y*dt;
dy = A(2,1)*x*dt + A(2,2)*y*dt;
dd = sqrt(dx.*dx + dy.*dy) + 0.01;  % Raw speeds, nonzero 
dy = dy./dd; dx = dx./dd; % Normalized unit arrows
quiver(x,y,dx,dy,ral);
title("Problem 3f");

% 3(g)  $\dot x=4x-3y$, $\dot y=8x-6y$.
A=[4 -3; 8 -6]; A, eig(A)
dx = A(1,1)*x*dt + A(1,2)*y*dt;
dy = A(2,1)*x*dt + A(2,2)*y*dt;
dd = sqrt(dx.*dx + dy.*dy) + 0.01;  % Raw speeds, nonzero 
dy = dy./dd; dx = dx./dd; % Normalized unit arrows
quiver(x,y,dx,dy,ral);
title("Problem 3g");

% 3(h)  $\dot x=y$, $\dot y=-x-2y$.
A=[0 1; -1 -2]; A, [L,V]=eig(A)
dx = A(1,1)*x*dt + A(1,2)*y*dt;
dy = A(2,1)*x*dt + A(2,2)*y*dt;
dd = sqrt(dx.*dx + dy.*dy) + 0.01;  % Raw speeds, nonzero 
dy = dy./dd; dx = dx./dd; % Normalized unit arrows
quiver(x,y,dx,dy,ral);
title("Problem 3h");
   
    






%Problem 4 (Ex.6.1.[8,9,10,11], p.199)
% Plot the flow fields in phase space for the following $2\times2$
% dynamical systems:
% (a) (Van der Pol's oscillator) $\dot x=y$, $\dot y=-x+y(1-x^2)$.
[x,y]=meshgrid(-4:4,-4:4);  % 2-D grid, integers, 9x9 
dx = y; dy = x+y.*(1-x.*x);  % 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);
title("Problem 4a = Ex.6.1.8 = Van der Pol");

% (b) (Dipole fixed point) $\dot x=2xy$, $\dot y=y^2-x^2$.
[x,y]=meshgrid(-4:4,-4:4);  % 2-D grid, integers, 9x9 
dx = 2*x.*y; dy = y.*y-x.*x;  % 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);
title("Problem 4b = Ex.6.1.9 = Dipole fixed point");

% (c) (Two-eyed monster) $\dot x=y+y^2$,
%  $\dot y=-\frac12 x+\frac15 y - xy +\frac65 y^2$.
[x,y]=meshgrid(-4:4,-4:4);  % 2-D grid, integers, 9x9 
dx = y+y.*y; dy = 0.5*x +0.2*y - x.*y + 1.2*y.*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);
title("Problem 4c = Ex.6.1.10 = Two-eyed monster");

% (d) (Parrot) $\dot x=y+y^2$, $\dot y=-x+\frac15 y -xy +\frac65 y^2$.
[x,y]=meshgrid(-4:4,-4:4);  % 2-D grid, integers, 9x9 
dx = y+y.*y; dy = -x +0.2*y - x.*y + 1.2*y.*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);
title("Problem 4d = Ex.6.1.11 = Parrot");



%Problem 5. (Ex.6.3.10, p.201)\
%    Consider the system $\dot x=xy$, $\dot y=x^2-y$, for which
%    linearization is inconclusive.
%
%    (a) Show that the linearization predicts that the origin is a non-isolated fixed point.
%
%    (b) Show that the origin is in fact an isolated fixed point.
%
%   (c) Classify the origin as attracting, repelling, saddle, or
%    something else, using a sketch of the flow field along the nullclines.
%
%    (d) Plot a computer-generated flow field phase portrait to confirm
%    your answer to part (c).

%%%% Flow field with nullclines
% 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.*y; dy = x.*x-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);
hold on;
ndxx=-4:0.1:4; ndxy=zeros(size(ndxx));  % Nullcline for dx
plot(ndxx,ndxy,"-.");
ndyx=-2:0.1:2; ndyy=ndyx.^2; % Nullcline for dy
plot(ndyx,ndyy,"--");
title("Problem 5d = Ex.6.3.10");
%% Some trajectories starting near (0,0)
% Rewrite the Xdot function for use in rk4a()
Xdot = @(V) [V(1).*V(2), V(1).^2-V(2)];
h=0.01;   % Use same time step for same accuracy
X = rk4a([ 0.1, 0.1], Xdot,  500,h); plot(X(:,1),X(:,2))
X = rk4a([-0.1, 0.1], Xdot,  500,h); plot(X(:,1),X(:,2))
X = rk4a([-0.1,-0.1], Xdot,  500,h); plot(X(:,1),X(:,2))
X = rk4a([ 0.1,-0.1], Xdot,  500,h); plot(X(:,1),X(:,2))
X = rk4a([ 0.4, 0.01], Xdot,  500,h); plot(X(:,1),X(:,2))
X = rk4a([-0.4, 0.01], Xdot,  500,h); plot(X(:,1),X(:,2))


%Problem 6  (Ex.6.4.[1,2,3], p.203)
% 6(a) $\dot x = x(3-x-y)$, $\dot y=y(2-x-y)$.
J=@(x,y)[3-2*x-y, -x ; -y , 2-x-2*y]
A=J(0,0); eig(A)
A=J(0,2); eig(A)
A=J(3,0); eig(A)
%  flow fields and nullclines
p=0:0.2:4; [x,y]=meshgrid(p,p);
dx=x.*(3-x-y); dy=y.*(2-x-y);
dd = sqrt(dx.*dx + dy.*dy) + 0.01;  % Raw speeds, nonzero
dy = dy./dd; dx = dx./dd; % Normalized unit arrows
clf; quiver(x,y,dx,dy,0.4)
hold on
p=0:0.1:3; plot(p,3-p);
p=0:0.1:2; plot(p,2-p);
title("6(a) Flow field and nullclines.")

% 6(b) $\dot x = x(3-2x-y)$, $\dot y=y(2-x-y)$.
J=@(x,y)[3-4*x-y, -x ; -y , 2-x-2*y]
A=J(0,0); eig(A)
A=J(0,2); eig(A)
A=J(3/2,0); eig(A)
A=J(1,1); eig(A)
%  flow fields and nullclines
p=0:0.2:4; [x,y]=meshgrid(p,p);
dx=x.*(3-2*x-y); dy=y.*(2-x-y);
dd = sqrt(dx.*dx + dy.*dy) + 0.01;  % Raw speeds, nonzero
dy = dy./dd; dx = dx./dd; % Normalized unit arrows
clf; quiver(x,y,dx,dy,0.4)
hold on
p=0:0.1:1.5; plot(p,3-2*p);
p=0:0.1:2; plot(p,2-p);
title("6(b) Flow field and nullclines.")


%6(c) $\dot x = x(3-2x-2y)$, $\dot y=y(2-x-y)$.
J=@(x,y)[3-4*x-2*y, -x ; -y , 2-x-2*y]
A=J(0,0); eig(A)
A=J(0,2); eig(A)
A=J(3/2,0); eig(A)
%  flow fields and nullclines
p=0:0.2:4; [x,y]=meshgrid(p,p);
dx=x.*(3-2*x-2*y); dy=y.*(2-x-y);
dd = sqrt(dx.*dx + dy.*dy) + 0.01;  % Raw speeds, nonzero
dy = dy./dd; dx = dx./dd; % Normalized unit arrows
clf; quiver(x,y,dx,dy,0.4)
hold on
p=0:0.1:1.5; plot(p,3/2-p);
p=0:0.1:2; plot(p,2-p);
title("6(c) Flow field and nullclines.")


%Problem 7 (Ex.6.5.2, p.207)
%  Consider the system $\ddot x=x-x^2$.
%  Convert to a first-order system:
% dx=y; dy=x-x^2=x(1-x).
%%%% (b) 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(-2:0.5:4,-4:0.5:4);  % 2-D grid, fine spacing 
dx = y; dy = x-x.*x;  % 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);
%%%% (c) Homoclinic orbit
hold on;
xh=0:0.01:3/2;
yh= sqrt(xh.^2-(2/3)*xh.^3);
plot(xh,yh); plot(xh,-yh);
title("Problem 7b,c = Ex.6.5.2");


%Problem 8 (Ex.6.5.4, p.207)
% Sketch the phase portrait for the system $\ddot x=ax-x^2$ for
%  $a<0$, $a=0$, and $a>0$.
% Convert to a first-order system:
% \dot x=y;\qquad \dot y=ax-x^2=x(a-x).
% Initialize the grid in phase space where arrows will
% indicate the flow X'=(x',y') at position X=(x,y):
[x,y]=meshgrid(-3:0.5:3,-2:0.5:2);  % 2-D grid, fine spacing 

a=-2;
dx = y; dy = a*x-x.*x;  % 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);
title("Problem 8 = Ex.6.5.4, a= -2");

a=0;
dx = y; dy = a*x-x.*x;  % 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);
title("Problem 8 = Ex.6.5.4, a=0");

a=2;
dx = y; dy = a*x-x.*x;  % 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);
title("Problem 8 = Ex.6.5.4, a=2");

%Problem 9 (Ex.6.6.2, p.212)
% System $\dot x=y$, $\dot y=x\cos y$ is reversible,
% Initialize the grid in phase space where arrows will
% indicate the flow X'=(x',y') at position X=(x,y):
[x,y]=meshgrid(-6:0.5:6,-4:0.5:4);  % 2-D grid, out past pi in y 
dx = y; dy = x.*cos(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);
title("Problem 9 = Ex.6.6.2");