function [temps,Y] = RUNGEKUTTA(t0,y0,T,N,f) 
% [temps,Y] = RUNGEKUTTA(t0,y0,T,N,f) 
h = T/N ; temps = [t0]; Y = [y0]; 
for k = 1:N 
  taux0 = temps(k); taux1 = temps(k)+h/2;  
  taux2 = temps(k)+h/2; taux3 = temps(k)+h; 
  y0    = f(taux0,Y(k)); 
  y1    = f(taux1,Y(k)+y0*h/2); 
  y2    = f(taux2,Y(k)+y1*h/2);  
  y3    = f(taux3,Y(k)+y2*h);  
  y     = Y(k) + h*(y0 +2*y1+2*y2 + y3)/6;  
  temps = [temps,taux3]; Y = [Y,y];  
end    
end