%% Example for system of ODEs

clearvars; clf; format short g

%% You first have to download vectorfield.m and put this file in the folder where your mlx files are.
%% Population of rabbits and foxes
% Let y_1(t) denote the population of rabbits, let y_2(t) denote the population of foxes at time t.
% We have for the rates of change y_1'(t),\; y_2'(t)
% Rabbits (prey)
%      y_1'(t) = R_1\cdot y_1  where R_1 = 2-y_2
% Foxes (predator)
%       y_2'(t)=R_2\cdot y  where R_2 = -1 + y_1
% We use the vector valued function \vec y(t) = \left[ \matrix{ y_1(t) \cr y_2(t) }
% \right] with derivative \vec y'(t) = \left[ \matrix{ y_1'(t) \cr y_2'(t) }
% \right]
% We can write the system of ODEs as
%       \vec y'(t) = \vec f(t,\vec y(t))
% Here we have
%        \vec f(t, \left[\matrix{y_1 \cr y_2}\right]) = \left[
% \matrix{ (2-y_2)y_1 \cr (-1+y_1)y_2 }
% \right]
% In Matlab we can define this function as

f = @(t,y) [ (2-y(2))*y(1) ; (-1+y(1))*y(2) ];

% We can evaluate e.g. \vec s = \vec f\Bigl(0,\left[\matrix{3\cr 4}\right]\Bigr) using

s = f(0,[3;4])

%% Initial value problem
% At the initial time t_0=0 we have initial populations y_1(0)=2,\; y_2(0)=2 for rabbits and foxes, i.e.
%       \vec y(t_0) = \vec y^{(0)} = \left[
% \matrix{2 \cr 2} \right]
%% Euler method
% We want to approximate \vec y(\frac12) using one step of the Euler method with h=\frac12

h = 1/2;

% At the initial time t_0 we have the value \vec y^{(0)}:

t0 = 0;
y0 = [2;2];

% We now evaluate \vec s = \vec f(t_0,\vec y^{(0)})

s = f(t0,y0)

% At the time t_1 = t_0+h=\frac12 we obtain the approximation \vec y^{(1)} = \vec y^{(0)} + h\cdot \vec s

y1 = y0 + h*s

% and obtain \vec y^{(1)}= \left[\matrix{2\cr 3}\right]. This means y_1(\frac12)\approx 2,\; y_2(\frac12)\approx 3.
%% Improved Euler method
% We now perform one step of the improved Euler method with h=\frac12
%       \vec s^{(1)} = \vec f(t_0,\vec y^{(0)})

s1 = f(t0,y0)

%       \vec y^E = \vec y^{(0)} + h\vec s^{(1)}

yE = y0+h*s

%       \vec s^{(2)} = \vec f(t_0+h,\vec y^E )

s2 = f(t0+h,yE)

%      \vec y^{(1)} = \vec y^{(0)} + h \frac12 \bigl(
% \vec s^{(s)} + \vec s^{(2)}
% \bigr)

y1 = y0 + h*(s1+s2)/2

%% Using ode45 to solve the initial value problem
% We want to approximate the solution \vec y(t) for t\in [t_0,T] with t_0=0,\;T=9:

[ts,ys] = ode45(f,[0,9],[2;2]);
disp([ts,ys]) % shows t, y1(t), y2(t) in each row (scroll to see all rows)

% Here ts is a vector of t-values, starting with t_0 and ending with T.
% ys is an array with two columns: the first column contains the approximations for y_1(t), the second column contains the approximations for y_2(t).
% The last row of ys contains the approximations for the final time T:

ys(end,:)         % last row of array ys, contains values for time T

% This means that y_1(9)\approx 1.8996, \; y_2(9)\approx1.6074.
% The first column ys(:,1) contains the population of rabbits. We can plot this using plot(ts,ys(:,1))

plot(ts,ys(:,1))
xlabel('t'); title('population y_1(t) of rabbits')

% The second column ys(:,2) contains the population of foxes.
% We can plot both functions y_1(t) and y_2(t) together using plot(ts,ys)

plot(ts,ys); legend('y_1','y_2'); xlabel('t')
title('y_1(t) (rabbits) and y_2(t) (foxes)')

% We start at a maximum of rabbits, then the foxes increase to a maximum, then the rabbits decrease to a minumum, then the foxes decrease to a minimum, then the rabbits increase to a maximum, etc.
% Note that the solution is periodic: after t\approx 4.6 everything repeats
%% Solution \vec y(t) in the "phase plane"
% The solution \vec y(t) = \left[ \matrix{y_1(t) \cr y_2(t)}\right] describes a point moving around in the y_1,y_2 plane, the so-called phase plane.
% The velocity vector of this moving point at time t is given by the vector \vec y'(t) = \vec f(t,\vec y(t)).
% In our example  the function \vec f(t,\vec y)=\left[
% \matrix{ (2-y_2)y_1 \cr (-1+y_1)y_2 }
% \right] does not depend on t, i.e., we have an autonomous ODE.
% We can plot the vector \vec f(t,\vec y) at a grid of points in the y_1,y_2 plane, this gives a "vector field":

vectorfield(f,0:.5:3,0:.5:4.25); hold on
xlabel('y_1 (rabbits)'); ylabel('y_2 (foxes)')

%
% Recall that ys(:,1) contains the rabbit population y_1(t) and ys(:,2) contains the fox population y_2(t).
% We can plot plot the points \bigl(
% y_1(t) , y_2(t)
% \bigr) in the phase plane using plot(ys(:,1),ys(:,2))

[ts,ys] = ode45(f,[0 9],[2;1]);
plot(ys(:,1),ys(:,2),'b')
hold off

% The periodic solution \vec y(t) corresponds to a closed curve in the phase plane.