%% Euler Method - Systems
% One can solve systems of ODEs using Euler's Method. Consider the predator-prey system from the text:
% \frac{dy}{dt} = 𝑦(1-𝛼𝑦)-\frac{𝑦𝑧}{1+𝛽𝑦}\\
% \frac{d𝑧}{𝑑𝑡}=-𝑧+\frac{𝑦𝑧}{1+𝛽𝑦}.
% We need to extend Euler's method to a vector form. First we modify the code as u_{i+1} = u_i + hf(u_i,v_i,t_i)\\
%   v_{i+1} = v_i + hg(u_i,v_i,t_i), where the system to be solved is \frac{du}{dt} =f(u,v,t),\\ \frac{dv}{dt} =g(u,v,t).

clear
u0 = [1; 0.01];
n=2000;
tspan = [0,80];      
[t,u] = euler2D(@predprey,tspan,u0,n);

%

y = u(:,1);  z = u(:,2);  
plot(t,y,t,z)
xlabel('t'), ylabel('u(t)'), title('Predator-prey solution')  
legend('prey','predator')

%
% This can be done more compactly using vecotrs. See the below code for explicit_euler.

[x,y] = explicit_euler(@predprey,tspan,u0,n); 
plot(x,y);
xlabel('t'), ylabel('u(t)'), title('Predator-prey solution')  
legend('prey','predator')
axis([tspan,0,9])

%
% One can compare this to the built-in function.

u0 = [1; 0.01];
t = linspace(0,80,2001)';       
[t,u] = ode45(@predprey,t,u0);

y = u(:,1);  z = u(:,2);  
plot(t,y,t,z)
xlabel('t'), ylabel('u(t)'), title('Predator-prey solution')  
legend('prey','predator')  

%
% Another example is a simple system using an anonymous function. Consider \frac{dx}{dt} = -0.2x+y\\
% \frac{dy}{dt} = -x-0.3y.

f = @(t,u) [-.2*u(1)+u(2);-u(1)-.3*u(2)];
tspan=[0,10];
u0=[1;0];
n=100;
[x,y] = explicit_euler(f,tspan,u0,n); 
plot(x,y);
xlabel('t'), ylabel('u(t)'), title('Solution Plots x(t) vs y(t)')  
legend('x(t)','y(t)')  

%
% Another way to visualize the solution is using a phase plane plot of y vs x.
%

plot(y(1,:),y(2,:));
xlabel('x(t)'), ylabel('y(t)'), title('Phase Plane') 
axis equal

%
% The functions used are below.

function dudt = predprey(t,u)
    alpha = 0.1;  beta = 0.25;
    y = u(1);  z = u(2);
    s = (y*z) / (1+beta*y);     % appears in both equations
    dudt = [ y*(1-alpha*y) - s;  -z + s ];  % RHS of system
end

%

function [t,U] = euler2D(dudt,tspan,u0,n)

a = tspan(1);  b = tspan(2);
h = (b-a)/n;
t = a + (0:n)'*h;
u = zeros(n+1,1);
v = zeros(n+1,1);
u(1) = u0(1);
v(1) = u0(2);
for i = 1:n
  du = dudt(t(i),[u(i),v(i)]);
  
  u(i+1) = u(i) + h*du(1,1);
  v(i+1) = v(i) + h*du(2,1);
end
U = [u,v];
end

function [x, y] = explicit_euler( f, xRange, y_initial, n )
x(1) = xRange(1);
h = ( xRange(2) - xRange(1) ) /n; 
[rows, cols] = size(y_initial);
y=zeros(rows,n+1);
y(:,1) = y_initial(:);
for k = 1 : n
    x(k + 1) = x(k) + h; 
    y(:,k+1) = y(:,k) + h * f( x(k), y(:,k) );
end
end