function [temps,X,Y] = PROIEPREDATEUR(param,epsilon,t0,x0,y0,T,N) 
a = param(1); b = param(2);
c = param(3); d = param(4);      
deff('[f]=f(x,y)','f = x*(a-epsilon*x-b*y)');
deff('[g]=g(x,y)','g = y*(-c+d*x)');     
h = T/N ; temps = [t0]; X = [x0]; Y = [y0]; 
for k=1:N 
    x0 = f(X(k),Y(k)); y0 = g(X(k),Y(k));  
    x1 = f(X(k) + x0*h/2,Y(k) + y0*h/2);  
    y1 = g(X(k) + x0*h/2,Y(k) + y0*h/2);  
    x2 = f(X(k) + x1*h/2,Y(k) + y1*h/2);  
    y2 = g(X(k) + x1*h/2,Y(k) + y1*h/2);  
    x3 = f(X(k) + x2*h,Y(k) + y2*h); 
    y3 = g(X(k) + x2*h,Y(k) + y2*h);
    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
endfunction 
