function [temps,X,Y] = PROIEPREDATEUR(param,epsilon,t0,x0,y0,T,N) 
% [temps,X,Y] = PROIEPREDATEUR(param,epsilon,t0,x0,y0,T,N) 
a = param(1); b = param(2);
c = param(3); d = param(4); 
f = inline('x*(u1-u0*x-u2*y)','x','y','u0','u1','u2'); 
g = inline('y*(-v1+v2*x)','x','y','v1','v2');      
h = T/N ; temps = [t0]; X = [x0]; Y = [y0]; 
for k=1:N 
    x0 = f(X(k),Y(k),epsilon,a,b); y0 = g(X(k),Y(k),c,d);  
    x1 = f(X(k) + x0*h/2,Y(k) + y0*h/2,epsilon,a,b);  
    y1 = g(X(k) + x0*h/2,Y(k) + y0*h/2,c,d);  
    x2 = f(X(k) + x1*h/2,Y(k) + y1*h/2,epsilon,a,b);  
    y2 = g(X(k) + x1*h/2,Y(k) + y1*h/2,c,d);  
    x3 = f(X(k) + x2*h,Y(k) + y2*h,epsilon,a,b); 
    y3 = g(X(k) + x2*h,Y(k) + y2*h,c,d);
    x  = X(k) + h*(x0+2*x1+2*x2+x3)/6; 
    y  = Y(k) + h*(y0+2*y1+2*y2+y3)/6;
    temps = [temps,temps(k)+h]; X = [X,x]; Y = [Y,y]; 
end
end