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

Horizon=100; // time Horizon

r0=0.05;rho=1/(1+r0);
// discount rate r0=5%

R=1/rho; // limit between sustainable and unsustainable
R=1.04; // unsustainable
R=1.1; // sustainable
growth=[1.04 1/rho 1.1];

for i=1:3
  R=growth(i);

// Construction of b and f through dynamic programming
b=[];
b(Horizon+1)=1;
for t=Horizon:-1:1
    b(t)=R*b(t+1)/(1+ R*b(t+1));
    end;

f=[];
f(Horizon+1)=0;
for t=Horizon:-1:1
    f(t)=(f(t+1)-log(rho*R))/(1+R*b(t+1));
    end;

// Optimal catches and stocks
opt=[];
hopt=[];
P0=10; // initial biomass
Popt(1)=P0;
for t=1:Horizon
     hopt(t)=min(Popt(t),max(0,b(t)*Popt(t)+f(t)));
    Popt(t+1)=max(0,R*(Popt(t)-hopt(t)));
    end

// Graphics
xset("window",10+i);xbasc();
plot2d2(0:(Horizon)',[[hopt; Popt($)] Popt],...
		rect=[0,0,Horizon*1.1,max(Popt)*1.2]);
// drawing diamonds, crosses, etc. to identify the curves
abcisse=linspace(1,Horizon,20); 
plot2d(abcisse,[ hopt(abcisse) Popt(abcisse) ],style=-[3,4]);
legends(['Catch trajectory';'Biomass trajectory'],-[3,4],'ur')
xtitle('Trajectories under linear growth with...
	   R='+string(R),'year (t)','biomass B(t)')

end

// \end{verbatim}