% 2016-09-10

% Logistic equation: 
% P' = r P (1 - P/K), P(0) = Po
% r>0, K>0 constants

clearvars; 
close all;

% Parameters

r = 1;
K = 1;
Po = 0.1;

% Time definitions

Tmax = 20;
dt = 0.01;
t = 0:dt:Tmax;

% Vector definitions

P = zeros(1,length(t));

% Initial conditions

P(1) = Po;

% Numerical computations

for j=1:length(t)-1
    k1p = r*P(j)*(1-P(j)/K);
    ap = P(j)+k1p*dt;
    k2p = r*ap*(1-ap/K);
    P(j+1) = P(j)+(k1p+k2p)*dt/2;
end

figure
hold on
plot(t,P,'-b','linewidth',2);
plot([0 Tmax],[0 0],'--','Color',[.7 .7 .7],'linewidth',1);
plot([0 Tmax],[K K],'--','Color',[.7 .7 .7],'linewidth',1);
axis([0 Tmax -0.1 1.1]);
set(gca,'fontsize',24);
xlabel('t  [a.u.]');
ylabel('P  [a.u.]');










