function [temps,Y] = HEUN(t0,y0,T,N,f) 
h = T/N ; temps = [t0]; Y = [y0]; 
for k = 1:N 
  taux = temps(k) + h; 
  yaux = Y(k) + h*f(temps(k),Y(k));      
  y    = Y(k) + h*(f(temps(k),Y(k))+f(taux,yaux))/2; 
  temps  = [temps,temps(k)+h]; Y  = [Y,y];
end    
endfunction   
