getf('plantI2.sci');
getf('plantII.sci');

T=10;
ateb=0.9;
r=1.2;
power = 0.5;

state=[0:(xm/6):(1.2*xp)];
control=0:0.05:1;

offspring=list();
discount=cumprod([1 ateb*ones(1,T-1)]);
for l=1:prod(size(control))
offspring(l)= (1-control(l))*(r*(state')^{power})*discount;
// offspring(l) is a matrix
end

final_cost_zero=zeros(state');



controlled_dynamics=dyn_plant_control;

[transition_matrix]=discretisation(state,control,controlled_dynamics);


instant_cost=offspring;
final_cost=final_cost_zero;

[value,feedback]=Bell_stoch(transition_matrix,...
instant_cost,final_cost);
// solves the dynamic programming equation

initial_state=grand(1,1,'uin',2,prod(size(state)));;
// corresponds to original state  >= state(2)
[z]=trajopt(transition_matrix,feedback,instant_cost,...
                                 final_cost,initial_state);
// computation of optimal trajectories (indexes)
zz=list();
zz(1)=state(z(1));
zz(2)=control(z(2));
zz(3)=z(3);
// optimal trajectories 

xset("window",1);xbasc();plot2d2(1:prod(size(zz(1))),zz(1));xtitle("taille")
xset("window",2);xbasc();plot2d2(1:prod(size(zz(2))),zz(2));xtitle("control")
xset("window",3);xbasc();plot2d2(1:prod(size(zz(3))),zz(3));xtitle("fitness")
// drawing of optimal trajectories