% fonction calculant la solution approche du systeme de Lotka-Volterra
%  s' = as-bsr
%  r' = -cr +dsr
% par le methode de Lotka Volterra de pas h
% valant r0 s0 en 0
% sur un intervalle [0,T]




function M=LVRK(s0,r0,a,b,c,d,T,h)


%le maillage de l'intervalle de temps
t=0:h:T;

%le nombre de noeuds du maillage.
Z= size(t);
z=Z(2);


% calcul des valeurs approchee par runge kutta
r(1)=r0;
s(1)=s0
for i=2:z
    k11=a*s(i-1)-b*s(i-1)*r(i-1);
    k12=-c*r(i-1)+d*s(i-1)*r(i-1);
    
    k21=a*(s(i-1)+h*k11/2)-b*(s(i-1)+h*k11/2)*(r(i-1)+h*k12/2);
    k22=-c*(r(i-1)+h*k12/2)+d*(s(i-1)+h*k11/2)*(r(i-1)+h*k12/2);

    k31=a*(s(i-1)+h*k21/2)-b*(s(i-1)+h*k21/2)*(r(i-1)+h*k22/2);
    k32=-c*(r(i-1)+h*k22/2)+d*(s(i-1)+h*k21/2)*(r(i-1)+h*k22/2);

    k41=a*(s(i-1)+h*k31/2)-b*(s(i-1)+h*k31/2)*(r(i-1)+h*k32/2);
    k42=-c*(r(i-1)+h*k32/2)+d*(s(i-1)+h*k31/2)*(r(i-1)+h*k32/2);


    s(i)=s(i-1)+h*(k11+2*k21+2*k31+k41)/6;
    r(i)=r(i-1)+h*(k12+2*k22+2*k32+k42)/6;
end

%On met les resultat dans une unique matrice M qui sortira les valeurs obtenues 
M(:,1)=t;
M(:,2)=s;
M(:,3)=r; 
