//\begin{verbatim}
// exec stab_BH.sce ; clear

////////////////////////////////////////////////////
// Antartic Blue whales data
////////////////////////////////////////////////////

// Biological parameters and functions

R_B=1.05; 
// Intrinsic rate of growth 5% per year
k = 400000; 
// Carrying capacity (whales)
b = (R_B-1)./k;

function [y]=Beverton(N)
// Beverton-Holt population dynamics
  y=(R_B*N)./(1 + b*N) 
endfunction

function [h]=Sust_yield(N)
// sustainable yield
  h=Beverton(N)-N
endfunction

function [y]=Beverton_e(t,N)
  he=Sust_yield(Ne);
  y=max(0,Beverton(N) - he);
endfunction


// Economic parameters

c_e=600000; // cost in $ per whale catcher year
p=7000; // price in $ per blue whale unit
q=0.0016; // catchability per whale catcher year



////////////////////////////////////////////////////
// Equilibria and yields
////////////////////////////////////////////////////

NMSE=  (sqrt(R_B) - 1)/b ;
// Maximum sustainable equilibrium MSE
NPPE= ( sqrt( R_B * (1 + (b*c_e/(p*q)) ) ) - 1 ) / b ;
// Private property equilibrium PPE
NCPE=c_e/(p*q) ;
// Common property equilibrium CPE

hMSE= Sust_yield(NMSE) ;
hPPE= Sust_yield(NPPE) ;
hCPE= Sust_yield(NCPE) ;


////////////////////////////////////////////////////
// Trajectories simulations 
////////////////////////////////////////////////////

// Only the last line counts. Change it.
Ne=NCPE     // unstable case
Ne=NPPE     // stable case 

he=Beverton(Ne)-Ne 

Horizon=500;  time=1:Horizon;
// Time horizon

epsil=20000;
// Perturbation level around the equilibrium 

xset("window",2); xbasc(2); 	
plot2d2(time,ones(1,Horizon)*NMSE,rect=[1,0,Horizon,k]);
Time=linspace(1,Horizon,20);
plot2d2(Time,ones(Time)*NMSE,style=-4);
legends("MSE",-4,'lr')
// Plot of the MSE

for (i=1:10)    // trajectories loop
  N0=Ne + 2*(rand(1)-0.5)*epsil;   
  // perturbation of initial state  
	Nt=ode("discrete",N0,0,time,Beverton_e); 
	// computation of the trajectory
	plot2d2(time,Nt,rect=[1,0,Horizon,k]);
	// plot of the trajectory
end

xtitle('Abundances trajectories','years',...
	   'abundances (whales)')

//\end{verbatim}