% Solution of the logistic growth with a threshold equation
% Numerical method: modified Euler (Runge-Kutta, order 2)

clearvars;
close all;

% ODE
% X' = -r (1-X/T) ( 1 - X/K) X + I;

r = 1;
T = 0.25;
K = 1;

I = 0.1;

logthr=@(x) -r*x.*(1-x/T).*(1-x/K);

% Computation of the fixed-point (assumes at most 3 fixed-points, brute force)

vv=-100:0.001:100;

fpcnt = 0;
if logthr(vv(1))+I > 0
    for j=1:length(vv)
        if logthr(vv(j))+I<0 && fpcnt == 0
            Vfp(1) = (vv(j)+vv(j-1))/2;
            fpcnt = 1;
        elseif logthr(vv(j))+I>0 && fpcnt == 1
            Vfp(2) = (vv(j)+vv(j-1))/2;
            fpcnt = 2;
        elseif logthr(vv(j))+I<0 && fpcnt == 2
            Vfp(3) = (vv(j)+vv(j-1))/2;
            fpcnt = 3;
        end
    end
end

% Computation of the eigenvalues (Note: they can be computed analytically)

dv = 0.001;
Eig = zeros(1,length(Vfp));
for k=1:length(Vfp)
    Eig(k) = (logthr(Vfp(k)+dv)-logthr(Vfp(k)-dv))/(2*dv);
end



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

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

x(1) = 0.26;

for j=1:length(t)-1
    k1x = logthr(x(j))+I;
    ax = x(j)+k1x*dt;
    k2x = logthr(ax)+I;
    x(j+1)=x(j)+(k1x+k2x)*dt/2;
end

xx = -1:0.01:2;

figure(1)
hFig = figure(1);
set(hFig, 'Position', [40 400 1000 500]); 
subplot(1,2,1)
hold on
plot(t,x,'-b','linewidth',2);
plot([0 Tmax],[T T],'--');
plot([0 Tmax],[0 0]','--');
plot([0 Tmax],[K K],'--');
axis([0 Tmax -0.1 1.1]);
set(gca,'fontsize',20);
xlabel('t');
ylabel('X');
subplot(1,2,2)
hold on
plot([xx(1)-1 xx(end)+1],[0 0],'--','Color',[.6 .6 .6]);
plot(xx,logthr(xx)+I,'linewidth',2);
axis([-0.1 1.1 -0.5 0.5]);
set(gca,'fontsize',20);
xlabel('X');
ylabel('F');

MVN = 0;
if MVN == 1
    xx = -1:0.01:2;
    figure(101)
    k = 1;
    hFig = figure(101);
    set(hFig, 'Position', [20 100 700 500])
    for j=1:floor(length(t)/k)-1
        subplot(1,2,1)
        plot(t,x,'-b','linewidth',2);
        hold on
        plot(t(k*j),x(k*j),'or','linewidth',2);
        plot([0 Tmax],[T T],'--');
        plot([0 Tmax],[0 0]','--');
        plot([0 Tmax],[K K],'--');
        axis([0 Tmax -0.1 1.1]);
        set(gca,'fontsize',20);
        xlabel('t');
        ylabel('X');
        pause(0.0001);
        hold off  
        subplot(1,2,2)
        plot([xx(1)-1 xx(end)+1],[0 0],'--','Color',[.6 .6 .6]);
        hold on
        plot(xx,logthr(xx)+I,'linewidth',2);
        plot(x(k*j),0,'or','linewidth',2);
        axis([-0.1 1.1 -0.5 0.5]);
        set(gca,'fontsize',20);
        xlabel('X');
        pause(0.0001);
        hold off   
    end
end

SMPFIG = 1;
if SMPFIG == 1
 

    % Sample figures

    close(1)

    figure
    hold on
    plot(t,x,'-b','linewidth',2);
    axis([0 Tmax -0.3 1.1]);
    set(gca,'fontsize',20);
    xlabel('t');
    ylabel('V');

    figure
    hold on
    plot([xx(1)-1 xx(end)+1],[0 0],'--','Color',[.6 .6 .6]);
    plot(xx,logthr(xx)+I,'linewidth',2);
    for k=1:length(Vfp)
        if Eig(k)<0
             plot(Vfp(k),0,'ok','linewidth',2);
             plot(Vfp(k),0,'xk','linewidth',2);
        elseif Eig(k)>0       
            plot(Vfp(k),0,'ok','linewidth',2);        
         end
    end
    axis([-0.3 1.1 -0.5 0.5]);
    set(gca,'fontsize',20);
    xlabel('V');
    ylabel('dV / dt');
    
   
end

