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

R_whale=1.05 ; 
// per capita productivity
R=R_whale ; 
K_whale = 400000;  
// carrying capacity  (BWH, blue whale unit)
K=K_whale ; 

// LINEAR DYNAMICS
function [y]=linear(B), y=R*B , endfunction;

Horizon= 10; 
years=1:Horizon;
yearss=1:(Horizon+1);

Binit= K/2 ; 
// initial condition
trajectory_whale=zeros(yearss) ;
// vector will contain the trajectory B(1),...,B(Horizon+1)
catch_whale=zeros(years) ;
// vector will contain the catches h(1),...,h(Horizon)

trajectory_whale(1)=Binit ;
// initialization of vector B(1),...,B(Horizon+1)
for t=years
  catch_whale(t)=R^{Horizon-t+1}/ ...
	  sum(R.^{Horizon-years(t:$)+1})*...
	  trajectory_whale(t);
  trajectory_whale(t+1)=...
	  linear(trajectory_whale(t)-catch_whale(t)); 
end

// Graphic display

xset("window",20+1); xbasc(20+1); 
plot2d2(yearss,[trajectory_whale ; [catch_whale 0]]' );
xtitle('Trajectories under linear growth with...
	   R='+string(R_whale),'year (t)',...
	   'biomass B(t) in blue whale unit')
plot2d2(yearss,[trajectory_whale ; [catch_whale 0]]',...
		style=-[3,4]);
legends(['Biomass trajectory';'Catch trajectory'],-[3,4],'ur')

// \end{verbatim}