%Euler's Method for solving logistic equation
% dN/dt = N(a-bN)
% N(t + dt) = N(t) + dt*f(N(t)) : Forward Euler

a = 2; %params for logistic
b = 1; %ditto
f = @(N) (N.*(a-b*N)); %function handle to RHS of ODE
N_sat = a/b; %saturation level
N0 = 0.1; %initial value e.g. initial population < saturation
NN0 = 5; %initial value e.g. initial population > saturation
t0 = 0; %initial time
t_final = 5; %final time

M = 100; %# of steps

dt = (t_final-t0)/M; %timestep
t = [0:dt:t_final]; %partition of time iterval

N = zeros(1,M+1); %initialize N
NN = zeros(1,M+1); %initialize N
N_sat_vec = N_sat*ones(1,M+1);

N(1) = N0; %initial value for N 
NN(1) = NN0; %initial value for NN

for i = 1:M
    N(i+1) = N(i) + dt*f(N(i)); %update for N, forward Euler
    NN(i+1) = NN(i) + dt*f(NN(i)); %update for NN, forward Euler
    figure(1);
    plot(t(1:i+1),NN(1:i+1),'k-x',t,N_sat_vec,'r-',t(1:i+1),N(1:i+1),'b-x')
    legend('N0 > a/b','saturation a/b','N0 < a/b')
    title('Numerical Approximation via Euler')
    
end

N_true = @(t) (a/b)./( 1 + ( (a-b*N0)/(b*N0) )*exp(-a*t)  );
NN_true = @(t) (a/b)./( 1 + ( (a-b*NN0)/(b*NN0) )*exp(-a*t)  );

figure(2)
plot(t,NN,'kx',t,N_sat_vec,'r-',t,N,'bx')
hold on;
plot(t,NN_true(t),'k-',t,N_true(t),'b-')
legend('approx (N0>a/b)','saturation','approx (N0<a/b)','actual','actual')
hold off
title('Actual vs Approximation')
    







