function [t,y] = Heun(t0,y0,T,N,f); 

% [t,y] = Heun(t0,y0,T,N,f); 
% Génère la méthode de Heun (formule des trapèzes) pour une EDO : y'=f(t,y)
% sur [t0,t0+T] x R. t0 est l'instant initial, y0 la valeur initiale T la
% longueur du segment d'etude, T/N le pas, f la fonction déclarée en ligne.
% En sortie, on récupère le temps t échantillonné et le vecteur
% correspondant y(t).

h = T/N; 
t = [t0]; 
y = [y0]; 
for k=1:N 
  A1 = f(t(k),y(k)); 
  A2 = f(t(k)+h,y(k)+h*A1); 
  yfinal = y(k) + (h/2)*(A1+A2);  
  t = [t,t(k)+h]; 
  y = [y,yfinal]; 
end  