% ### ODErkEX1.m ###     09.19.14
% numerically solve the Logistic equation
% P'(t)= k*P(1-P/L)

% Program calculates in four ways: 1&2. solve explicitly using Euler and RK4, 
% 3. solve via ode45 (via external function call) and 4. actual solution

clear; figure(1); clf

% ************************
% User Inputs

k= 1;   % intrinsic growth rate (const.)
L= 5;   % carrying capacity (const.)
Pinit= L/15;   % initial condition at tI(1)
tI= [0 10];  % time boundaries
stepsize= 0.05;     % for RK4

% ************************
% Runge-Kutta 4 (and Euler's method)
m=1;  % counter
for j=tI(1):stepsize:tI(2)
    if j == tI(1)
        P(m)= Pinit;    Peuler(m)= Pinit;
    else
        P0= k*P(m-1)*(1-P(m-1)/L);  % deriv. at y=y0 (last meas.)
        P1= k*(P(m-1)+(stepsize/2)*P0)*(1-(P(m-1)+(stepsize/2)*P0)/L);
        P2= k*(P(m-1)+(stepsize/2)*P1)*(1-(P(m-1)+(stepsize/2)*P1)/L);
        P3= k*(P(m-1)+(stepsize)*P2)*(1-(P(m-1)+(stepsize)*P2)/L);
        P(m)= P(m-1) + (stepsize/6)*(P0+ 2*P1+ 2*P2+ P3);   % RK4 solution
        Peuler(m)= P(m-1) + stepsize* P0;      % also store away Euler's method value
    end
    t(m)= j;  % keep track of t for plotting
    m= m+ 1;  % increment counter 
end

% visualize
plot(t,Peuler,'g+'); grid on;   hold on
plot(t,P,'o--'); 
xlabel('t'); ylabel('P(t)'); 

% ************************
% can also solve via Matlab's builtin ode45, but need to use an external
% function to define the ODE
[tM,PM] = ode45(@(t,P) logistic(t,P,k,L),tI, Pinit);
plot(tM,PM,'kx','LineWidth',2);

% ************************
% also plot analytic solution
A= (L-Pinit)/Pinit;
sol= L./(1+A*exp(-k*t));
plot(t,sol, 'r-')
legend('Euler','RK4','ode45 solution','Exact solution','Location','SouthEast')
title('Solution to logisitc equation using various numerical methods')