//\begin{verbatim}
// exec dyn_pop.sce

lines(0);

r = 0.9; K= 10; R = 1+r ; MVP=K/5;
// Population dynamics parameters  

k=R*K/(R-1);
function [y]=Logistic(t,B)
  y=max(0, (R*B).*(1-B/k) )
endfunction
            
function [y]=Ricker(t,B)
  y=B.*exp(r*(1-B/K))
  endfunction
       
  b=(R-1)/K ; 
function [y]=Beverton(t,B)
  y=(R*B)./(1 + b*B) 
endfunction

function [y]=Depensation(t,B)
//  y=max(0,B+(Beverton(t,B)-B).*(B-MVP)/MVP)
    y=max(0,B+(Beverton(t,B)-B).*(B-MVP))
  endfunction
// Dynamics

xset("window",0);xbasc(0); 
B=linspace(0,2*K,1000);
plot2d(B,[B' Ricker(0,B)' Logistic(0,B)' Beverton(0,B)'...
		  Depensation(0,B)']);
// drawing diamonds, crosses, etc. to identify the curves
B=linspace(0,2*K,30);
plot2d(B,[B' Ricker(0,B)' Logistic(0,B)' Beverton(0,B)'...
		  Depensation(0,B)'],-[1,2,3,4,5]);
legends(["Identity";"Ricker";"Logistic";"Beverton";...
		 "Depensation"],-[1,2,3,4,5],'ul');
xtitle('Biomass dynamics','Biomass B(t)','Biomass B(t+1)')
// Comparaison of the shapes of population dynamics
      

T=10; time=0:T; 
// Time horizon

N_simu=50;
// N_simu=30;
// Number of simulations

xset("window",1:4); xbasc(1:4); 
// opening windows

for i=1:N_simu
// simulation loop

	B_0=rand(1)*1.5*K;
	// random initial conditions
	
	y_Ricker=ode("discrete",B_0,0,time,Ricker);
	y_Logistic=ode("discrete",B_0,0,time,Logistic);
	y_BH=ode("discrete",B_0,0,time,Beverton);
	y_D=ode("discrete",B_0,0,time,Depensation);
	// Computation of trajectories starting from B0 
		// along distinct population dynamics 
	
	xset("window",1); 
	plot2d(time,[y_Ricker'],rect=[0,0,T,2*K]);
	xtitle('Ricker Trajectories','time (t)',...
	'biomas B(t)')
		
	xset("window",2); 
	plot2d(time,y_Logistic,rect=[0,0,T,2*K]);
	 xtitle('Logistic Trajectories','time (t)',...
	 'Biomass B(t)')
		 
	xset("window",3); 
	plot2d(time,[y_BH'],rect=[0,0,T,2*K]);
	xtitle('Beverton-Holt Trajectories','time (t)',....
	'Biomass B(t)')
		
	xset("window",4); 
	plot2d(time,[y_D'],rect=[0,0,T,2*K]);
	xtitle('Depensation Trajectories','time (t)',...
	'Biomass B(t)')
       
end
// end simulation loop

//\end{verbatim}