// Exemple algorithme de Howard 

function M= Mtrans(n)
  M=rand(n,n); 
  s=sum(M,'c'); 
  M= M./(s*ones(1,n)); 
endfunction

// Howard Pb de temps d'arret 

function y=phi(x) 
  y = - x^2 
endfunction

// les états 
n=5; 
M=Mtrans(n);
xs=(1:n)' 
r=0.02;
// controle initial: le max est réalisé par phi 
u = 2*ones(xs);
v = 0*ones(xs);
vnew=v;

Mc = [M,zeros(n,1);zeros(1,n),1];
Ms = [0*M,ones(n,1);zeros(1,n),1];
u  =[u;1];
v  =[v;0];
b  =[phi(xs);0];
cs =b;
cc =0*b;

while %t 
  // u donné: résoudre 
  // v = (1/(1+r))phi la ou u=2 
  // v = (1/(1+r))Mv  la ou u=1 
  ic = find(u== 2);
  is   = find(u== 1);
  M1=Mc; 
  M1(is,:)=Ms(is,:); 
  vnew= ((1+r)*eye(M1) -M1)\ b 
  // On a calculé un nouveau v 
  write(%io(2),[mini(vnew-v),maxi(vnew-v)])
  v = vnew;
  // On calcule un nouveau u 
  [m,u]=min(Ms*v+cs,Mc*v+cc);
end