//\begin{verbatim}

// -----------------------------------------
// DYNAMIC PROGRAMMING ALGORITHM
// -----------------------------------------

VV=zeros(MM,T); // value functions in a matrix 
// The final value function is zero.
UUopt=(NN+1)*ones(MM,T); // optimal controls in a matrix

for t=(T-1):(-1):1 // backward induction
  loc=[];
  // will contain the vector to be maximized
  loc_bottom=mini(VV(:,t+1));
  // The value attributed to the value function VV(:,t)
  // when a control is not admissible.
  //
  for uu=1:NN // columns 1 to NN selected
	loc=[loc, Admissible(:,uu) .* ( ...
	   discount^t * instant_gain(:,uu) + ...
	  VV( Future_Integer(:,uu) , t+1 ) ) +...
	+ (1-Admissible(:,uu)) .* (discount^t * bottom + loc_bottom) ] ;
  end
// When the control uu is admissible, 
// loc is the usual DP expression.
// When the control uu is not admissible, 
// loc is the  DP expression
// with both terms at the lowest values.
	 //
  loc=[loc, VV(:,t+1) + discount^t *0] ;
// Adding an extra control/column which provides zero 
// instantaneous gain and does not modify the state:
// retire option.
//
  [lhs,rhs]= maxi(loc,"c") ; // DP equation
  VV(:,t)=lhs;
  UUopt(:,t)=rhs;
//
UUopt(Stop_Integers,t)=(NN+1)*ones(Stop_Integers) ;
	// retire option   
end


// -----------------------------------------
// OPTIMAL TRAJECTORIES
// -----------------------------------------

xx=zeros(T,NN);
zz=zeros(T,1);
uu=(NN+1)*ones(T,1);
vv=0;

xx(1,:)=zeros(1,NN);
// initial profile

xset("window",1) ; 
xset("colormap",hotcolormap(64));
// table of colors
xbasc();plot2d2(-0.5+0:(NN+2),-[0 xx(1,:) 0 0],...
				rect=[0,-HH-0.5,NN+1,0.5])
Matplot1(rich_color,[0,-HH-0.5,NN+1,0.5]);

t=1 ; last_control=0 ;
//
while last_control<NN+1 do 
  zz(t)=profile2integer(xx(t,:)) ;
  uu(t)=UUopt(zz(t),t) ;
  xx(t+1,:)=xx(t,:);
  if uu(t)<NN+1 then
  xx(t+1,uu(t))=xx(t,uu(t))+1;
  vv=vv+discount^t*richness(xx(t+1,uu(t)),uu(t)) ;
  end
//  halt()
  xbasc();plot2d2(-0.5+0:(NN+2),-[0 xx(t+1,:) 0 0],...
				  rect=[0,-HH-0.5,NN+1,0.5]); xpause(400000);
  
  last_control=uu(t);
t=t+1 ;
end
xtitle("Optimal extraction profile for discount rate "...
	   +string(discount))


function []=display_profile(time)
  // Displays graphics of mine profiles
xset("window",time) ; 
xset("colormap",hotcolormap(64));
xbasc();
plot2d2(-0.5+0:(NN+2),-[0 xx(time,:) 0 0],...
		rect=[0,-HH-0.5,NN+1,0.5]);
// Just to to set the frame
Matplot1(rich_color,[0,-HH-0.5,NN+1,0.5]);
// The value of the mine blocks in color. 
// Will be in background.
plot2d2(-0.5+0:(NN+2),-[0 xx(time,:) 0 0],...
		rect=[0,-HH-0.5,NN+1,0.5]);
xtitle("Optimal extraction profile at ultimate time t="...
	   +string(time) +" for discount rate "+string(discount))
endfunction 

time=t;
display_profile(time)

time=1;
display_profile(time)

time=ceil(T/3);
display_profile(time)

//\end{verbatim}