function [T,Y]=outbreak(xo,param_ini,param_fin,tfin)
%
% Demonstrate the insect outbreak
%
%  Usage: [T,Y]=outbreak(x0,[k_ini,r_ini],[k_fin,r_fin],tfin); 
%
%  Example:
%      
%clear
%global r_ini r_fin k_ini k_fin tfin aa;
%k_ini=1;k_fin=5;r_ini=0.5;r_fin=1;
%k_ini=15;k_fin=15;r_ini=0.6;r_fin=0.025;
%k_ini=15;k_fin=15;r_ini=0.2;r_fin=0.8;
%%k_ini=15;k_fin=15;r_ini=0.530;r_fin=0.555;
%tfin=5000;
%[T,Y]=outbreak(1,[k_ini,r_ini],[k_fin,r_fin],tfin);
%

global r_ini r_fin k_ini k_fin tfin aa;
format compact;

options= odeset('Reltol',1.e-5);
%options= odeset('Reltol',1.e-3,'Outputfincn','odephas3','Refine',4);
lw=2;
fs=16;
fs2=19;

 k_ini=param_ini(1);
 k_fin=param_fin(1);
 r_ini=param_ini(2);
 r_fin=param_fin(2);
 IC=[xo];
 tspan=[0 tfin];
 %aa=0.001;
 aa=5/tfin;
 
 [T,Y]=ode45(@outbreak_eqns,tspan,IC,options);
 
 R=(r_ini+(r_fin-r_ini)/2)+(r_fin-r_ini)*tanh(aa*(T-tfin/2))/2;
 K=(k_ini+(k_fin-k_ini)/2)+(k_fin-k_ini)*tanh(aa*(T-tfin/2))/2;
 
 % ploting the orbit 
 figure(1); clf
 set(gca,'FontSize',fs)
 subplot(2,1,1)
 plot(T,Y,'-o','LineWidth',lw)
 xlabel('t','FontSize',fs2); ylabel('x','FontSize',fs2)
 title('Insect Outbreak','FontSize',fs2)
 set(gca,'FontSize',fs)
 subplot(2,2,3)
 plot(T,R,'-','LineWidth',lw)
 axis tight
 hold on
 plot(T,R*0+0.5373,'--r','LineWidth',lw)
 plot(T,R*0+0.2628,'--r','LineWidth',lw)
 xlabel('t','FontSize',fs2); ylabel('r','FontSize',fs2)
 set(gca,'FontSize',fs)
 subplot(2,2,4)
 plot(T,K,'-','LineWidth',lw)
 xlabel('t','FontSize',fs2); ylabel('k','FontSize',fs2)
 set(gca,'FontSize',fs)
 hold off
 
 figure(2); clf
 set(gca,'FontSize',fs)
 xx=[1.025:0.01:3,3:21];
 subplot(1,1,1)
 rr=2*xx.^3./((1+xx.^2).^2);
 kk=2*xx.^3./(xx.^2-1);
 plot(kk,rr,'--r','LineWidth',lw)
 hold on
 plot(K,R,'-bo','LineWidth',lw)
 plot(K(1),R(1),'*r','LineWidth',lw)
 axis tight;
 axis([0 40 0 0.8]);
 xlabel('k','FontSize',fs2); ylabel('r','FontSize',fs2)
 title('Parameter Diagram','FontSize',fs2)
 set(gca,'FontSize',fs)


 figure(3);clf
 x=[0:0.25*2:25];
 r=[0:0.01*2:1];
 k=[0.25:0.125*2:18];
 [RR,KK,XX]=meshgrid(r,k,x);
 RHS=RR.*(1-XX./KK)-XX./(1+XX.^2);
 plotisosurface(RR,KK,XX,RHS,0,'red');
 view([-50 35])
 hold on
 mm=max(x)/10000;
 plot3(R,K,Y+mm,'bo-')
 plot3(R,K,Y-mm,'bo-')
 hold off

 xlabel('r','FontSize',fs2); 
 ylabel('k','FontSize',fs2)
 zlabel('x','FontSize',fs2)
 title('Orbit','FontSize',fs2)
 set(gca,'FontSize',fs)
 box on
 grid on
 axis tight
 return;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Differential equations defining the Lorenz system
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function F=outbreak_eqns(t,vec)
global r_ini r_fin k_ini k_fin tfin aa;

 %r=r_ini+(r_fin-r_ini)*t/tfin;
 %k=k_ini+(k_fin-k_ini)*t/tfin;
 
 r=(r_ini+(r_fin-r_ini)/2)+(r_fin-r_ini)*tanh(aa*(t-tfin/2))/2;
 k=(k_ini+(k_fin-k_ini)/2)+(k_fin-k_ini)*tanh(aa*(t-tfin/2))/2;

 x=vec(1);

 F = r*x*(1-x/k)-x^2/(1+x^2);

 return;

