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

R_tuna=2.25; // discrete-time intrinsic growth
R=R_tuna;
K_tuna = 250000;  // carrying capacity in metric tons
K = K_tuna;
c_tuna=2500; 
// unit cost of effort in dollars per standard fishing day
p_tuna=600; // market price in dollars per metric ton 

// BEVERTON-HOLT DYNAMICS 
R_BH = R_tuna ; 
b_BH = (R_BH-1) / K ;

// SUSTAINABLE YIELD FUNCTION
function [SY]=sust_yield(B)
  SY=B-( B ./ ( R_BH - b_BH *B ) ) ; 
endfunction


B_MSE = ( R_BH - sqrt(R_BH) ) / b_BH ;
// maximum sustainable equilibrium

MSY=sust_yield(B_MSE) ; 
// maximum sustainable yield

xset("window",1); xbasc();
abcisse=linspace(0,K,100);
plot2d(abcisse,sust_yield(abcisse),rect=[0,0,K,1.1*MSY]);
H_MSY=linspace(0,MSY,20);
plot2d(B_MSE*ones(H_MSY),H_MSY,style=-6);
xtitle("Sustainable yield for the ...
	   Beverton-Holt model (tuna)",...
	   "biomass (metric tons)","catches (metric tons)");

// private property equilibrium
cost=c_tuna;
price=p_tuna;
B_PPE = ...
	( R_BH - sqrt(R_BH - (b_BH * cost / price) ) ) / b_BH ;

//\end{verbatim}