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

% [t,y] = RungeKutta(t0,y0,T,N,f); 
% Génère la méthode de Runge-Kutta classique (formule de Simpson) 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).

% on declare d'abord le pas :  
h = T/N ; 
% on initialise le vecteur des temps : 
t = [t0];  
% on initialise le vecteur de y 
% (y(t0)=y0 figure ici la condition initiale) : 
y = [y0]; 
% On redige maintenant la boucle de calcul suivant 
% les instructions decrites ci-dessus en (5.37) : 
for k=1:N
  % on declare les points intermediaires entre t0+(k-1)*h et t0+k*h : 
  t1 = t(k);  
  t2 = t1 + h/2 ; 
  t3 = t1 + h/2 ; 
  t4 = t1 + h ; 
  % on calcule les quatre valeurs auxiliaires A1, A2, A3, A4 : 
  A1 = f(t1, y(k)); 
  A2 = f(t2, y(k) + h*A1/2); 
  A3 = f(t3, y(k) + h*A2/2); 
  A4 = f(t3, y(k) + h*A3) ;
  % on calcule y(t0+(k+1)*h) comme indique :  
  yfinal  = y(k) + (h/6) *(A1 + 2*A2 + 2*A3 + A4) ; 
  % on complete les vecteurs des temps et des y :
  t  = [t,t4] ;  
  y  = [y,yfinal] ;
end 