// On definit la fonction
//-------------------------

k1 = 0.04;
k2 = 300000;
k3 = 1000;

// On definit les donnees 
//-----------------------
h    = 0.01
xx   = [0:h:1];
yy   = [0:h:0.5];
epsi = 1.E-5;

// On trace les isovaleurs de f et g pour avoir un aperçu 
// de leur comportement.
//--------------------------------------------------------

function res = F(x,y);
    res = -k1*x + k3*y*(1.-x-y); 
endfunction

function res = G(x,y);
    res = +k1*x - k3*y*(1-x-y) - k2*y*y;  
endfunction

//disp("On trace les isovaleurs de la fonction f(x,y)");
//zz1 = feval(xx,yy,F);
//xset('window',1)    // ouvre la fenêtre 1
//contour(xx,yy,zz1,10,zlev=0.01)

//disp("On trace les isovaleurs de la fonction g(x,y)");
//zz2 = feval(xx,yy,G);
//xset('window',2)    // ouvre la fenêtre 2
//contour(xx,yy,zz2,10,zlev=0.01)

// On resout le syst. diff. en appliquant un solveur implicite
//------------------------------------------------------------------


// On definit des fonctions
//--------------------------

function res = F3(Yp)
res = [-k1*Yp(1)+k3*Yp(2)*Yp(3); k1*Yp(1)-k3*Yp(2)*Yp(3)-k2*Yp(2)*Yp(2); k2*Yp(2)*Yp(2)];    
endfunction

function res = dF3(Yp)
    res = [-k1, k3*Yp(3), k3*Yp(2); k1,-k3*Yp(3)-2*k2*Yp(2),-k3*Yp(2);0,2*k2*Yp(2),0];
endfunction

function resu = MonNewton(dt,Y0)
  Y1=Y0;
  erreur=1.;
  while (erreur > epsi) 
      A = eye(3,3) - dt*dF3(Y1);
      Y2 = Y1 - A\(Y1-dt*F3(Y1)-Y0);
      erreur=norm(Y1-Y2)
      Y1=Y2;
  end
  resu=Y1
endfunction    


// On applique la méthode 
//------------------------

Y0=[1;0;0]
disp("Y0 vaut",Y0);

T  = 40.
tn = 0;
dt = 1.E-2;
sot= 100;

N=int(T/(sot*dt));
solu = ones(3,N);
tt   = ones(1,N);

solu(:,1)=Y0;
tt(1) = 0.;

it=2;
nn=0;

while (tn < T)
    Y1 = MonNewton(dt,Y0);
    //Y1 = Y0+dt*F3(Y0); 
    tn = tn+dt;
    if modulo(nn,sot) == 0 then
     solu(:,it) = Y1;
     tt(it)     = tn; 
     it = it+1;
     mass = Y1(1) + Y1(2) + Y1(3);
     disp("la masse vaut",mass);
    end
    nn=nn+1;
    Y0=Y1;
    
    

end

xset("line style",2); plot2d(tt,solu(1,:),style=5);
xset("line style",4); plot2d(tt,solu(2,:),style=3);
xset("line style",6); plot2d(tt,solu(3,:),style=1);
legends(["yA";"yB";"yC"],[[5;2],[3;4],[1;6]],opt="lr");


